I need to numerically perform an integral of the form

where F is a certain function, \mathbf{u}(t) is the solution of a simple ODE (it is a 1D Hamiltonian system) and \frac{\partial \mathbf{u}(t)}{ \partial \mathbf{u}(0)} is the Jacobian of \mathbf{u}(t) with respect to the initial condition \mathbf{u}(0). For that reason, I need a fast way of evaluating \mathbf{u}(t) and \frac{\partial \mathbf{u}(t)}{ \partial \mathbf{u}(0)} for many values of t.

The following code allows me to calculate these values, but it isnâ€™t efficient:

```
using DifferentialEquations, ForwardDiff
function f!(du,u,p,t)
du[1] = exp(-u[2])*(1-exp(-u[2]))
du[2] = 2*u[1]
end
function get_trajectory(u0,t_max)
problem = ODEProblem(f!,u0,(0.0,t_max))
solve(problem)
end
function get_jacobian(u0,t_max)
problem(u) = ODEProblem(f!,u,(0.0,t_max))
t->ForwardDiff.jacobian(u->solve(problem(u))(t),u0)
end
u = get_trajectory([1.0,0.0], 1.0)
u(0.5)
jac = get_jacobian([1.0,0.0], 1.0)
jac(0.5)
```

It works fine for calculating \mathbf{u}, as, when I call `get_trajectory`

, the differential equation is solved once and then each call of `u(t)`

is fast. On the other hand, each time I call `jac(t)`

the differential equation is solved once again, which makes the call slow. I was wondering if there is a way to solve the differential equation only once and get all the values I need.

On the DifferentialEquations.jl documentation there is a page on Sensitivity Analysis that appears to treat a similar question, but the Jacobian is calculated with respect to parameters. I also found this topic which may have a solution to my problem, but the links containing the references to the solution appear broke to me.