Hi

I’m still in the newbie category so I dare ask this question as a “First Steps” question.

I read this blog:

```
https://www.johndcook.com/blog/2020/02/08/arenstorf-orbit/
```

… about an equation that was used to construct orbits for the Apollo program. It looked interesting, and not too complex, so I thought I could implement that in Julia. So I tried, and failed. The program is “numerically unstable” Julia tells me.

Dear Julia community, can you at a glance see what I’m doing wrong? I’ve tried the DPRKN12, DPRKN6 and KahanLi8 solvers, but they seem to be failing in more or less the same way so it’s probably something other than the solver that is wrong, but I have no idea what.

```
using DifferentialEquations
using OrdinaryDiffEq
using Plots
function plot_arensdorf()
μ = 0.012277471
μ′ = 1 - μ
function arensdorf_orbit(du, u, z, p, t)
x, y = z
x′,y′ = u
D(z) = ((x + z )^2 + y^2)^(3/2)
D1 = D(μ)
D2 = D(-μ′)
x′′ = x + 2y′ + μ′*(x + μ)/D1 - μ*(x - μ′)/D2
y′′ = y + 2x′ + μ′*y/D1 - μ*y/D2
du[1] = x′′
du[2] = y′′
end
initial_positions = [0.994, 0]
initial_velocities = [0.0, -2.001585106]
tspan = (0.0, 5000.0)
prob = SecondOrderODEProblem(arensdorf_orbit, initial_velocities, initial_positions, tspan)
sol = solve(prob, DPRKN12(), dt=1/100);
plot(sol,vars=(1,2))
end
```

```
julia> plot_arensdorf()
┌ Warning: Instability detected. Aborting
└ @ DiffEqBase ~/.julia/packages/DiffEqBase/YIwj5/src/integrator_interface.jl:339
Error showing value of type Plots.Plot{Plots.GRBackend}:
ERROR: InexactError: trunc(Int64, -Inf)
Stacktrace:
[1] trunc at ./float.jl:702 [inlined]
[2] ceil at ./float.jl:365 [inlined]
[3] optimize_ticks_typed(::Float64, ::Float64, ::Bool, ::Array{Tuple{Float64,Float64},1}, ::Int64, ::Int64, ::Int64, ::Float64, ::Float64, ::Float64, ::Float64, ::Bool, ::Nothing) at /Users/rmz/.julia/packages/PlotUtils/6AioP/src/ticks.jl:183
[4] #optimize_ticks#36 at /Users/rmz/.julia/packages/PlotUtils/6AioP/src/ticks.jl:137 [inlined]
[5] (::PlotUtils.var"#kw##optimize_ticks")(::NamedTuple{(:k_min, :k_max),Tuple{Int64,Int64}}, ::typeof(optimize_ticks), ::Float64, ::Float64) at ./none:0
[6] optimal_ticks_and_labels(::Plots.Subplot{Plots.GRBackend}, ::Plots.Axis, ::Nothing) at /Users/rmz/.julia/packages/Plots/qZHsp/src/axes.jl:188
[7] optimal_ticks_and_labels at /Users/rmz/.julia/packages/Plots/qZHsp/src/axes.jl:156 [inlined]
[8] get_ticks(::Plots.Subplot{Plots.GRBackend}, ::Plots.Axis) at /Users/rmz/.julia/packages/Plots/qZHsp/src/axes.jl:266
[9] _update_min_padding!(::Plots.Subplot{Plots.GRBackend}) at /Users/rmz/.julia/packages/Plots/qZHsp/src/backends/gr.jl:902
[10] iterate at ./generator.jl:47 [inlined]
[11] _collect(::Array{AbstractLayout,2}, ::Base.Generator{Array{AbstractLayout,2},typeof(Plots._update_min_padding!)}, ::Base.EltypeUnknown, ::Base.HasShape{2}) at ./array.jl:635
[12] collect_similar(::Array{AbstractLayout,2}, ::Base.Generator{Array{AbstractLayout,2},typeof(Plots._update_min_padding!)}) at ./array.jl:564
[13] map(::Function, ::Array{AbstractLayout,2}) at ./abstractarray.jl:2073
[14] _update_min_padding!(::Plots.GridLayout) at /Users/rmz/.julia/packages/Plots/qZHsp/src/layouts.jl:310
[15] prepare_output(::Plots.Plot{Plots.GRBackend}) at /Users/rmz/.julia/packages/Plots/qZHsp/src/plot.jl:261
[16] display(::Plots.PlotsDisplay, ::Plots.Plot{Plots.GRBackend}) at /Users/rmz/.julia/packages/Plots/qZHsp/src/output.jl:143
[17] display(::Any) at ./multimedia.jl:323
[18] #invokelatest#1 at ./essentials.jl:709 [inlined]
[19] invokelatest at ./essentials.jl:708 [inlined]
[20] print_response(::IO, ::Any, ::Bool, ::Bool, ::Any) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.3/REPL/src/REPL.jl:156
[21] print_response(::REPL.AbstractREPL, ::Any, ::Bool, ::Bool) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.3/REPL/src/REPL.jl:141
[22] (::REPL.var"#do_respond#38"{Bool,REPL.var"#48#57"{REPL.LineEditREPL,REPL.REPLHistoryProvider},REPL.LineEditREPL,REPL.LineEdit.Prompt})(::Any, ::Any, ::Any) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.3/REPL/src/REPL.jl:719
[23] #invokelatest#1 at ./essentials.jl:709 [inlined]
[24] invokelatest at ./essentials.jl:708 [inlined]
[25] run_interface(::REPL.Terminals.TextTerminal, ::REPL.LineEdit.ModalInterface, ::REPL.LineEdit.MIState) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.3/REPL/src/LineEdit.jl:2306
[26] run_frontend(::REPL.LineEditREPL, ::REPL.REPLBackendRef) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.3/REPL/src/REPL.jl:1045
[27] run_repl(::REPL.AbstractREPL, ::Any) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.3/REPL/src/REPL.jl:201
[28] (::Base.var"#770#772"{Bool,Bool,Bool,Bool})(::Module) at ./client.jl:382
[29] #invokelatest#1 at ./essentials.jl:709 [inlined]
[30] invokelatest at ./essentials.jl:708 [inlined]
[31] run_main_repl(::Bool, ::Bool, ::Bool, ::Bool, ::Bool) at ./client.jl:366
[32] exec_options(::Base.JLOptions) at ./client.jl:304
[33] _start() at ./client.jl:460
julia>
```