Yeah this comes up from time to time and it seems legit to me. Here’s a previous related post: Costs/benefits of solving large ODE systems as systems of difference equations/in discrete time (and I think that links to another one before it)
In this case I agree with @wds15 it’s probably the scaling of the forward sensitivity problem that is messing things up.
When you implement the ODE solver in Stan with a fixed timestep you don’t have a sensitivity scaling problem. Similarly if you use the prototype adjoint sensitivity thing then the scaling should be better, but then it could still end up being slower if trapezoid is a more appropriate solver than the BDF or Adams solvers in CVODES.
There’s some solver functions @jtimonen wrote for @spinkney 's helpful stan functions repo here (look under functions/odes). The RK4 method there might be faster than the trapezoid method.