# How to construct the variational equation of an ODE without handwriting the Jacobi matrix

Hello! I am new to Julia. I want to solve the variational equation of an ODE without handwriting the Jacobi matrix. Here is an MME. Suppose that the ODE equations are the followings:

\dot{x}=y\\ \dot{y}=x-x^3-0.1y+0.1\cos(t)

Set f(x,y,t)=(y,x-x^3-0.1y+0.1\cos(t))^T. Given a initial condition (x_0,y_0) at t=0, suppose the solution is (x(t),y(t)) for t\in[0,2\pi]. Now I want to solve the following matrix ODE:

\dot{X}=d_{(x,y)}f(x(t),y(t),t)X

with initial condition X(0)=id. I can solve the above matrix ODE by using the following codes:

using OrdinaryDiffEq, StaticArrays

function vectorfield(u, p, t)
SA[u[2], u[1]-u[1]^3-0.1*u[2]+0.1*cos(t)] # vector field of the system
end

function timemap(u0)
tspan = (0.0, 2 * pi)
prob = ODEProblem(vectorfield, u0, tspan)
sol = solve(prob, Vern7())
function dv(X, p, t)
SA[0.0 0.1; 1-3*(sol(t)[1])^2 -0.1] * X # variational equation in matrix form
end
X0 = SA[1.0 0.0; 0.0 1.0] # initial condition of variational equation
probJ = ODEProblem(dv, X0, tspan)
solJ = solve(probJ, Vern7(), save_everystep=false)
(sol[end], solJ[end])
end


However, when the ODE has a large dimension, hand calculation of the Jacobi matrix is not a good choice. Can anyone show me how to solve the matrix ODE like the above example without hand calculation of the Jacobi matrix?

What you’re looking for seems to be automatic (or algorithmic) differentiation. How large is “large dimension”? Depending on the answer, there can be several packages for computing the Jacobian: see the list at https://juliadiff.org/