Hessian of function involving integration of ODE

Hi all,
I apologize if it’s a possibly stupid question, but I am trying to calculate the Hessian of a function that internally performs the integration of a differential equation. What I am currently doing is this:

using Zygote, DifferentialEquations
function der!(du, u, p, t)
    du[1] = -p[1] * u[1]

function cost_for_hessian(params)
    prob_uode_pred = ODEProblem{true}(der!, [10.0], (0, 10))
    solutions = solve(prob_uode_pred, Tsit5(), p=params, saveat=[10.0],  sensealg=ForwardDiffSensitivity())
    return sum(abs2.((solutions)))

sensitivity_analysis_hessian = hessian(p -> cost_for_hessian(p), [2.0])

When I try to execute the code (julia version 1.9.1, Zygote v0.6.61, DifferentialEquations v7.7.0), I get the following error: ERROR: Cannot determine ordering of Dual tags.

Could anyone tell me if what I am trying to do is feasible and how to do it?
Thank you in advance!

Use Forward over Reverse, not Zygote.hessian. These examples

make use of that. It’s done in Optimization.jl like this:

Or you can directly use Direct Adjoint Sensitivities of Differential Equations · SciMLSensitivity.jl

Thank you!