Hello!

I was working on the following problem where I’m integrating a second order differential equation. Initially I tried putting it into the 2nd Order Solver in DifferentialEquations.jl, where it did produce a plot set, but when comparing to colleagues they were not correct.

Later, I changed the 2nd Order ODE System to a First Order ODE System, which this one does produce the correct results. I’m able to work with the First Order System, but now I’m more curious as to why my implementation of the 2nd Order version is producing different results. For context, these are the Hill equations. I tried different integrators for the `solve()`

, and some did change the plots, but none reproduced the plots the first order system created.

I provided the code block beneath that should show the two different cases.

Any understanding on why the results are different would be appreciated!

```
using DifferentialEquations
using Plots
function main()
b = range(1, 3, length=6)
t0 = -100.
t1 = 150.
tspan = (t0, t1)
# system_ode(b, t0, tspan)
second_ode(b, t0, tspan)
end
function system_ode(b, t0, tspan)
function hill!(du, u, p, t)
f = -3/(u[1]^2 + u[2]^2)^1.5
du[1] = u[3]
du[2] = u[4]
du[3] = u[1]*(f + 3) + 2u[4]
du[4] = u[2]*f - 2u[3]
end
index = size(b, 1)
orbits = Array{Any}(nothing, index)
u0 = zeros(4)
for i in 1:index
u0[1] = b[i] - 8/(3b[i]^2 * (t0^2 + 4/9)^0.5)
u0[2] = -1.5*u0[1]*t0
u0[3] = 0
u0[4] = -1.5*u0[1]
prob = ODEProblem(hill!, u0, tspan, abstol=1e-8, reltol=1e-8)
sol = solve(prob)
orbits[i] = plot(sol[1,:], sol[2,:], label="b = "*string(b[i]), framestyle=:origin, xlim=(-2,2), ylim=(-2,2))
end
plot(orbits..., layout=index)
end
function second_ode(b, t0, tspan)
function hill!(ddu, du, u, p, t)
f = -3/(u[1]^2 + u[2]^2)^1.5
ddu[1] = u[1]*(f + 3) + 2du[2]
ddu[2] = u[2]*f - 2du[1]
end
index = size(b, 1)
orbits = Array{Any}(nothing, index)
u0 = zeros(2)
du0 = zeros(2)
for i in 1:index
u0[1] = b[i] - 8/(3b[i]^2 * (t0^2 + 4/9)^0.5)
u0[2] = -1.5*u0[1]*t0
du0[1] = 0
du0[2] = -1.5*u0[1]
prob = SecondOrderODEProblem(hill!, du0, u0, tspan, abstol=1e-8, reltol=1e-8)
sol = solve(prob)
orbits[i] = plot(sol[1,:], sol[2,:], label="b = "*string(b[i]), framestyle=:origin, xlim=(-2,2), ylim=(-2,2))
end
plot(orbits..., layout=index)
end
main()
```