@DavidKetcheson's answer hits the big points: you can always construct methods of high enough order using extrapolation, that's a very pessimistic bound and you can always do a whole lot better, all the good ones are derived by hand (with the help of some computer algebra tools), no lower bound is known, and the highest order methods are due to Feagin. Given some of the comments, I wanted to round out the answer with a discussion of the current state-of-the-art tableaus in the field.
If you want a compendium of RK tableaus, you can find one in this Julia code. Citations for the paper that they came from are in the docstrings for the tableau constructors. The developer documentation for DifferentialEquations.jl lists all of these tableaus as available for use, and you can see here that these are all tested using Travis and AppVeyor continuous integration suites to make sure that not only the order conditions are satisfied, but they they actually achieve the requested convergence (verification testing). From these, you can see that there are:
- 5 order 9 methods
- 6 order 10 methods
- 2 order 12 methods
- 1 order 14 method
(that I could find that were published). Again, all derived by hand.
The convergence tests show that some of the derivations were not carried out to high enough precision to work for more than 64-bit numbers (they are commented like this). So that's an interesting quirk to be aware of: at these high orders you usually only get coefficients that "to an error
x" satisfy the order conditions, but when using arbitrary precision arithmetic you can actually detect these bounds. So the precision to which you carry out the coefficients do matter, and you should choose it to cover the precision you wish to test (/ use, of course).
If you want a bunch of stability plots, you can just
plot(tableau) using the Plots.jl recipe. A good set of notes which has a lot of this written down can be found on Peter Stone's website (go below and click on say the order 10 schemes and you'll get a bunch of PDFs). When developing DifferentialEquations.jl, I created that set of tableaus to systematically go through them on test problems / look at the analytical indicators to see which ones should be included in the main library. I made some quick notes here. As you can see from the algorithms included in the main library, the ones that I found worthwhile were the Verner and Feagin methods. The Verner 9th order method is the highest order method with an interpolant matching the order too. That's something to recognize: the Feagin methods do not have a matching interpolant (though you can bootstrap Hermite, but that's really inefficient).
Since they are all implemented with very efficient implementations, you can play around with them yourself and see how much the different features actually matter. Here's a Jupyter notebook which shows the Feagin methods in use. Notice that the convergence plot really is going to
1e-48 error. High order methods are only more efficient than lower order methods when you really need a very very low tolerance. You can find some benchmarks that use some of them at DiffEqBenchmarks.jl, though when they are used it's usually the 9th order Verner method, and usually showing that the benchmark is not in the regime where this high of order is efficient.
So if you want to play around and work with some high order methods, RK-Opt is what I found is a good tool for deriving some (as @DavidKetcheson mentioned), and DifferentialEquations.jl has all of the published methods (I think?) implemented so that you can easily test/benchmark against them. However, unless you find an assumption that can be dropped, from my tests I haven't been able to find something that beats the Verner (orders 6-9) and Feagin (orders 10+) methods. YMMV though, and I'd love to see more research in this.