I have a system of ~10 ODEs with ~20 parameters and a loss function (e.g. MSE) that is evaluated on the trace of the ODE solution. What’s an efficient way to calculate the Hessian matrix of the loss w.r.t the parameters? Is there an autodiff package that works with ODE solvers that allows second derivatives? Would appreciate a link to an example if possible!