For medium sized stiff system biology ODE models I have recently done two benchmarks that suggest that `CVODE_BDF`

outperform the Julia solvers (in particular the multi-step bdf solvers).

Below is an with 25 ODEs and 39 parameters. The model is imported with SBMLImporter.jl and I used PEtab.jl to set the parameters to the values in the original publication. The model files can be found here. I ran the benchmarks using Julia 1.10.3

```
using DiffEqBase, OrdinaryDiffEq, Catalyst, SBMLImporter, Sundials, Plots, DiffEqDevTools,
TimerOutputs, LinearAlgebra, ModelingToolkit, BenchmarkTools, PEtab
# Read as ODESystem
prnbng, cb = load_SBML(joinpath(@__DIR__, "Bachmann_MSB2011", "model_Bachmann_MSB2011.xml"),
check_massaction=false)
osys = convert(ODESystem, prnbng.rn)
# Read as PEtabODEProblem and extract parameters for first simulations condition
petab_model = PEtabModel(joinpath(@__DIR__, "Bachmann_MSB2011", "Bachmann_MSB2011.yaml"), verbose=false)
petab_prob = PEtabODEProblem(petab_model, verbose=false)
_p = PEtab.get_ps(petab_prob.θ_nominalT, petab_prob)
_u0 = PEtab.get_u0(petab_prob.θ_nominalT, petab_prob)
tf = 100.0
tspan = (0.0, tf)
oprob = ODEProblem{true, SciMLBase.FullSpecialize}(osys, _p, tspan, _u0)
oprob_jac = ODEProblem{true, SciMLBase.FullSpecialize}(osys, _p, tspan, _u0, jac=true)
oprob_sparse = ODEProblem{true, SciMLBase.FullSpecialize}(osys, _p, tspan, _u0, jac=true, sparse=true)
sol = solve(oprob, CVODE_BDF(), abstol=1/10^14, reltol=1/10^14)
test_sol = TestSolution(sol)
abstols = 1.0 ./ 10.0 .^ (6:10)
reltols = 1.0 ./ 10.0 .^ (6:10)
setups = [
Dict(:alg=>CVODE_BDF()),
Dict(:alg=>TRBDF2()),
Dict(:alg=>QNDF()),
Dict(:alg=>FBDF()),
Dict(:alg=>KenCarp4()),
Dict(:alg=>Rodas4()),
Dict(:alg=>Rodas5P())
];
wp = WorkPrecisionSet(oprob, abstols, reltols, setups; error_estimate=:l2,
saveat=tf/10., appxsol=test_sol, maxiters=Int(1e9), numruns=100)
plot(wp)
title!("No Jacobian")
wp_jac = WorkPrecisionSet(oprob_jac, abstols, reltols, setups; error_estimate=:l2,
saveat=tf/10., appxsol=test_sol, maxiters=Int(1e9), numruns=100)
plot(wp_jac)
title!("With Jacobian")
setups[1] = Dict(:alg=>CVODE_BDF(linear_solver=:KLU))
wp_sparse = WorkPrecisionSet(oprob_sparse, abstols, reltols, setups; error_estimate=:l2,
saveat=tf/10., appxsol=test_sol, maxiters=Int(1e9), numruns=100)
plot(wp_sparse)
title!("Sparse Jacobian")
```

With output:

Is there any additional performance option I am missing for the Julia solvers?