In the figure below you can see that ODEInterfaceDiffEq.radau5 and RadauIIA5 show the same behavior; they are wrong for tol=1e-12 but correct at tol=1e-14. But at t~500 seconds they are more than 20 times slower than VCABM. I am not familiar with using BigFloat. I will figure it out and try using bigfloats to see whether there is an advantage. Is there a disadvantage when it comes to speed when using BigFloats? In Python I have seen a huge degradation in performance when an arbitrary precision library like mpmath is used.
Edit: Big float at tol=1e-12 just slows down the code by a factor of 2 and at low tolerance gives wrong results. The plot below is for tol=1e-6.
Edit2: Vern9 does not run. It says instability detected: aborting. Is there a way to not abort Vern9 when instability is detected?