Just to clarify here, I am pretty sure both are possible (I use AD on an ODE solver all the time), see https://docs.sciml.ai/DiffEqDocs/v5.3/analysis/sensitivity.html https://arxiv.org/abs/1812.01892 AD is quite useful here because you can get the sensitivity of non-analytic things easily. E.g. sensitivity of gravity assist periapsis wrt control. Also using AD on ODE solvers saves a lot of time because you don't need to write out the equations to propagate the Jacobian, they are implicitly created with the dual numbers. As something to watch out for when doing this, you need to be really careful about adaptive timestepping. A timestep which is valid for the original problem may not be valid for the problem with AD - the timestepping algorithm needs to take this into account.