>As for Sundials and FatODE, my understanding was that they used finite differences, forward mode, or differentiating the solver operations for at least some aspect of their sensitivity analysis.

No, their sensitivity analysis doesn't use finite differences. They perform the sensitivity analysis as described in their documentation. If methods for the Jacobian calculation are not provided, then they utilize finite differences on the Jacobian calculation of course, using the same routine as for the stuff solver. But with a given Jacobian function there's no finite differences or forward mode.

>This isn't true - when we do the reverse solve, we get all gradients using a _single_ solve going backwards in time. I'm now realizing that the misunderstanding might be caused by our Fig. 2, which shows that multiple backward solves are necessary when the loss depends on the state at multiple time points.

No. Of course it's a single solve going backwards in time. I see what my misread was, but I'm surprised you'd attempt to solve the equation backwards like that because it's known to not be stable. Without a reversible integrator (implicit Adams is not reversible), it's a well-known result that the method drifts from the true solution doing a backwards integration, so the values so z(t) needs to computed with forward passes in order to be correct. A good test equation for this is probably the Lornez equation with standard parameters over a time like [0,300]. The backwards pass will diverge and not necessarily be on the same butterfly wing as the forwards pass. The Julia code using CVODE is shown here ( https://gist.github.com/ChrisRackauckas/fef4ae7778320530d44b... ) and you see that starting from [1.0,1.0,1.0] the result of going backwards is then off by [-17.5445, -14.7706, 39.7985]. So there you go, CVODE's Adams method ends up on a different "wing" of the butterfly when integrated backwards, ending up not even close to the actual initial point (CVODE is the successor to LSODE, both by Alan Hindmarsh, but utilizes constant leading coefficient forms to reduce computations. So not exactly the same as the paper, but very close). Thus to ensure correctness, existing sensitivity analysis packages only get Jacobians of f using data from forward passes. The Neural ODEs may have had a small enough Lyopunov coefficient or a short enough integration that this wasn't an issue, but it is in general something to note. Of course, if you are only doing this on Hamiltonian systems...

>I agree that our statement about the stability of implicit over explicit was overly broad, thanks for pointing that out. Can you suggest a more accurate statement about the advantages of implicit over explicit methods?

Any statement on it is too broad to be useful. Runge-Kutta Chebyshev methods are explicit methods for stiff systems. Implicit Adams is an implicit method for non-stiff systems. And there's many more examples. The stiffness handling also depends on implementation details. Using functional iteration on a BDF method reduces the region of stability, which is why BDF needs to use Newton's method for solving the implicit equation in order to be applicable to stiff ODEs. It's best to just talk about the stability of individual methods and their implementation.