Yes, there aren't too many resources on this for some reason. For a very long time, the standard goto was "just use BDF methods". This mantra was set in stone for few historical reasons: for one Gear's code was the first widely available stiff solver, and for another the MATLAB suite didn't/doesn't include an implicit RK method. However, this heuristic isn't always correct, and I would say from testing it's usually wrong. Let me explain in detail.
BDF methods have a high fixed cost
Adaptive timestepping and adaptive order in BDF methods has a really high cost. Coefficients have to be recalculated, or values need to be interpolated to different times. There has been a lot of work going into making the current BDF codes do better here (there are two well-known "forms" for implementation trying to handle this issue differently), but it's actually just a very difficult software engineering problem. Because of this, actually most of the BDF codes are all descendants from Gear's original code: Gear's, vode, lsoda, Sundial's CVODE, even the DAE solvers DASKR and DASSL are descendants from the same code.
What this means is that, if you problem is "too small", the high fixed cost really matters and Implicit RK methods will do better.
High order BDF methods are very unstable for complex eigenvalues
BDF methods allow you to control the maximum order and make it more conservative for a reason: the higher order BDF methods cannot handle even "medium sized" complex eigenvalues very well. So in these cases, in order to be stable, they have to drop their order. This is the reason why the 6th order BDF method, while technically stable, is often ignored because even the 5th order already has issues here (and the 6th order is even less stable). Only up to 2nd order is A-stable, so it can always fallback to there, but then stepping is error dominated.
So when using BDF codes on non-trivial problems, you're not getting 5th order all of the time. Oscillations cause its order to drop.
BDF methods have a high starting cost
BDF methods which are autostarting have a high starting cost. They start with a small Euler step, then a small BDF-2 step, etc. For shorter time integrations, this has a very non-negligible effect. If you are stopping often, say for event handling, this will heavily impair the efficiency of the codes. This is one reason why you rarely if ever see multistep methods used in heavy event handling or delay equation situations (delay equations restart at each high enough order discontinuity, which is propagated at each $i\tau$ for every delay-time $\tau$)
BDF methods, being multistep methods, have the best scaling
But this is what BDF methods have going for them: they scale the best. This is because they have the least function calls. This erroneously has led to the idea that you should "always default to BDF methods" (which for the reason above is not always true for stability and efficiency reasons), but it does mean that, if the function call $f$ is expensive enough, then BDF methods will be best. Generally this means that for large enough and stiff enough PDE discretizations, BDF methods will be the most efficient. That said, this "large enough" is very problem dependent, and I have found PDE discretizations which do better with radau methods of up to 1,000's of spatial variables, so YMMV. Only way to know is to test.
What benchmarks are available?
Hairer's book and DiffEqBenchmarks (explained below) are probably the best in terms of easily available work-precision diagrams. Hairer's Solving Ordinary Differential Equations II has a bunch of work-precision diagrams on pages 154 and 155. The results on the problems he chose match what I stated above for the reasons I stated above: implicit RK will be more efficient if the problem is not "sufficiently large". Another interesting thing to note is that high order Rosenbrock methods come out as the most efficient in many of his tests (like Rodas) in the regime where error is higher, and the implicit RK radau5 is the most efficient at lower errors. But his tests problems are mostly not discretizations of large PDEs, so take the points above into account.
How do you test/benchmark?
I like to test this with DifferentialEquations.jl in Julia (disclaimer: I am one of the developers). This is in Julia. Programming language should really get a note here. Remember that as the cost of the function call increases, BDF methods fair better. In R/MATLAB/Python, the user's function is the only R/MATLAB/Python code that the optimized solvers have to actually use: if you're using SciPy or Sundials wrappers, it's all C/Fortran code except for the function you pass. This means that, in dynamic languages (which aren't Julia), BDF methods do better than they normally would because the function calls are very unoptimized.
So Julia is a good testing ground for the more optimal case since if the function call is type-stable, your $f$ is as fast/efficient as any C/Fortran function. What's nice with DifferentialEquations.jl is you can test between a ton of algorithms just by switching one line of code. For an ODEProblem prob, you can switch between the two with:
@time sol = solve(prob,CVODE_BDF())
@time sol = solve(prob,radau())
where the first uses Sundials' CVODE (BDF method), and the second uses Hairer's radau (implicit RK).
For any ODEProblem, you can use the benchmarking tools to see how the different algorithms scale for different adaptive tolerances. Some results are posted at DiffEqBenchmarks.jl. For example, on the ROBER problem (system of 3 stiff ODEs) you can see that Sundials actually is unstable and diverges with a high enough tolerance (whereas the other methods converge just fine), showing the note above about stability issues. On the Van Der Pol problem, you can see it's more of a wash. I have a lot more than I haven't posted, but probably won't get to it until I finish up some higher order Rosenbrock methods (Rodas is the Fortran version of those).
(Note: those stiff benchmarks need updating. For one, the text needs updating since for some reason ODE.jl's methods diverge...)
Extrpolation methods and RKC
There are also extrapolation methods like seulex which are made for stiff problems. These are "infinite adaptive order", but that only means they are asymtopically good when you're looking for very low error (i.e. looking to solve stiff problems at lower than 1e-10 or so, in which case you can probably use an explicit method though). However, in most cases they are not as efficient and should be avoided.
Runge-Kutta Chebyschev is an explicit method which also works on stiff problems that you should consider. I don't have it wrapped in DifferentialEquations.jl yet so I don't have any hard evidence myself of how it fairs.
Other methods to consider: specialized methods for stiff PDEs
Should probably make a note that high order Rosenbrock methods do really well, many times even better, for small-medium sized stiff problems when you can easily compute the Jacobian. However, for some PDEs, I believe advection-diffusion problems fall into this category, Rosenbrock can lose some orders of accuracy. Also, they need very accurate Jacobians in order to keep their accuracy. In Julia this is easy because the solvers come with symbolic and autodifferentiation which can be correct to machine epsilon. However, things like MATLAB's ode23s can have issues because these implementations use finite differencing. For BDF and implicit RK methods, errors in the Jacobian lead to slower convergence, while for Rosenbrock, since these are not implicit equations and are rather RK methods with Jacobian inversions in there, these just lose order of accuracy.
Other methods to look at are exponential RK methods, exponential time-differencing (ETD), implicit integration factor (IIF), and exponential Rosenbrock methods. The first three make use of the fact that, in many PDE discretizations,
$$ u_t = Au + f(u) $$
where $A$ is a linear operator, and this linearity allows you to exactly solve the $Au$ operator part. The downside is, the resulting method has to use $e^A$ which is dense even if $A$ was sparse, and so there's a lot of work in Krylov methods. Still, these are A-stable methods which avoid implicit equations and so they can be very efficient for PDE discretizations.
Exponential Rosenbrock methods don't need the $A$, and instead split by $Ju + g(u)$ where $J$ is the Jacobian, and $f = Ju + g$, but it's the same idea.
Still other methods: the latest creations
Methods which are fully implicit obviously do well for stiff equations. If the PDE is not large enough, you cannot effectively "parallelize in space" enough to make use of HPCs. Instead, you can create parallel-in-time discretizations which are fully implicit and thus good for stiff equations, yet parallel to make full use of hardware. XBraid is a solver that uses a technique like this, and the main methods are PFFAST and parareal methods. DifferentialEquations.jl is developing a neural net method which works similarly.
Again, this is optimal when you don't have a large enough spatial discretization to efficiently parallelize, but have the resources for parallel computation available.
Conclusion: take asymptotic considerations with a grain of salt
For a very long time in numerical diffeq history, there was the idea that extrapolation methods were the best because they can technically achieve any order. This means that, for a small enough $\Delta t$, they will be the most efficient (so, if you're looking for a very precise solution, yes!). However, almost no one actually works in the regime where these are efficient, which is why people generally don't use extrapolation methods anymore. This is a word of caution: take asymptotic considerations with a grain of salt
BDF methods are asymptotically the best, but in most cases you probably aren't working in that regime. But if the spatial discretization has enough points, BDF methods can efficiently parallelize in space (by parallelizing the linear solving) and will have the least function calls, and thus will do the best. But if your PDE discretization is not large enough, take a good look at implicit RK, Rosenbrock, exponential RK, etc. methods. Using a software suite like DifferentialEquations.jl that makes it easy to swap between the different methods can be really helpful for you to understand what method does best for your problem domain, since in many cases it cannot be known in advance.
[If you have any example problems you want added to the testing suite, please feel free to help get something in there. I want to keep a very comprehensive resource on this. I hope to have all of Hairer's benchmarks in runnable notebook forms "soon", and would like other problems which are representative of scientific fields. Any help is appreciated!]
