using DifferentialEquations, ForwardDiff

function lotka_volterra!(du, u, p, t)
  x, y = u
  α, β, δ, γ = p
  du[1] = dx = α*x - β*x*y
  du[2] = dy = -δ*y + γ*x*y
end

# Initial condition
u0 = [1.0, 1.0]

# Simulation interval and intermediary points
tspan = (0.0, 10.0)
tsteps = 0.0:0.1:10.0

# LV equation parameter. p = [α, β, δ, γ]
p = [1.5, 1.0, 3.0, 1.0]

# Setup the ODE problem, then solve
prob = ODEProblem(lotka_volterra!, u0, tspan, p)
sol = solve(prob, Tsit5(), saveat=tsteps)

function solvep(p)
    _prob = remake(prob,p=p)
    solve(_prob, Tsit5(), saveat=tsteps,abstol=1e-12,reltol=1e-12)[1,:]
end

using ForwardDiff
ForwardDiff.jacobian(solvep,p)


function jacf(u,p,t)
    x, y = u
    α, β, δ, γ = p
    [α-β*y -β*x
     γ*y   -δ+γ*x]
end

function lotka_volterra_sense!(du, u, p, t)
    x, y = u
    α, β, δ, γ = p
    du[1] = dx = α*x - β*x*y
    du[2] = dy = -δ*y + γ*x*y
    jac = jacf(u,p,t)
    du[3:4]  = jac * u[3:4]  + [x,0]
    du[5:6]  = jac * u[5:6]  + [-x*y,0]
    du[7:8]  = jac * u[7:8]  + [0,-y]
    du[9:10] = jac * u[9:10] + [0,x*y]
end
u0sense = [1.0, 1.0, zeros(8)...]

probsense = ODEProblem(lotka_volterra_sense!, u0sense, tspan, p)
solsense = solve(probsense, Tsit5(), saveat=tsteps, abstol=1e-12,reltol=1e-12)

jac1 = solsense[[3,5,7,9],:]'
jac2 = ForwardDiff.jacobian(solvep,p)

jac1 - jac2