I realize this is not a very good question due to the nature of the problem, but I have to try since I’ve spent hours debugging this. I am interested in solving a set of ODEs that can be written as follows:

```
using Elliptic.Jacobi
function dn_ab!(du,u,p,t)
du[1] .= +im*p[1].*u[1] .+ im.*u[2].*dn.(t,real(p[2]))
du[2] .= -im*p[1].*u[2] .+ im.*u[1].*dn.(t,real(p[2]))
end
# Initial Condition
function ab_dn_t0(x)
A = +χ .+ 0.5.*Ω*λ.*(x) .- π/4
B = -χ .+ 0.5.*Ω*λ.*(x).- π/4
return [2*im*sin(A); 2*cos(B)]
end
```

I asked about this before here: Solving same DiffEq repeatedly but with different initial conditions

The equation is solved multiple times and the solutions are recursed over and combined nonlinearly to generate a final solution, \psi, which is the only thing I have some physical intuition about. I will include a MWE showing how this \psi is generated (sorry about that, I did my best to simplify things, I just have no intuition at all on how to deal directly with the solutions of the ODE since they are just intermediate steps).

Using any algorithm except `Euler()`

to solve the equation, I get a garbage solution, something like (using `Tsit5()`

):

The correct solution, only obtainable using `Euler()`

.

Everything else used to generate this final solution \psi is correct, I have spent countless hours checking, and the result is perfectly correct when I use `Euler()`

.

My question are:

- Why does nothing work but
`Euler()`

- Is there a way to use a different algorithm? Some flag I missed in the documentation somewhere?

MWE (I did my best to make it as small as possible, I promise):

You only need to change the line near the end that defines `algo`

to see the drastic effect on the solution.

I’ve tried almost all the explicit non-stiff algorithms in the documentation and they all give various types of garbage results.

```
using Elliptic.Jacobi
using DifferentialEquations
using Plots
function solve!(λ, m, χ, Ω, Nₓ, Nₜ, x, t, algo)
N = length(λ)
ψₜ = Array{Complex{Float64}}(undef, Nₜ, Nₓ, N+1)
ψₜ[:,:,1] .= exp.(im*x'.*(1-m/2)).*dn.(t,m)
calc_rs(N, 1, ψₜ, λ, χ, Ω, Nₓ, Nₜ, x, t, algo)
return ψₜ[:,:,N+1]
end
function calc_rs(n, p, ψₜ, λ, χ, Ω, Nₓ, Nₜ, x, t, algo)
@info "Calculating Lax pair generating functions rₙₚ(x,t) and sₙₚ(x,t) for (n,p) = ($n,$p)"
# Base case
if n == 1
# Allocate empty arrays to hold lax pair generating functions
rf = Array{Complex{Float64}}(undef, Nₜ, Nₓ)
sf = Array{Complex{Float64}}(undef, Nₜ, Nₓ)
# Transpose x
x = x'
# Set up a funnction for the ODE
function dn_ab!(du,u,p,t)
du[1] = +im*p[1].*u[1] .+ im.*u[2].*dn.(t,real(p[2]))
du[2] = -im*p[1].*u[2] .+ im.*u[1].*dn.(t,real(p[2]))
end
# And the initial condition
function ab_dn_t0(x)
A = +χ[p] .+ 0.5.*Ω[p]*λ[p].*(x) .- π/4
B = -χ[p] .+ 0.5.*Ω[p]*λ[p].*(x) .- π/4
return [2*im*sin(A); 2*cos(B)]
end
# Prepare the ODE intnegrator
tspan = (0.0,abs(minimum(t)))
u0 = ab_dn_t0(x[1])
prob = ODEProblem(dn_ab!, u0, tspan, [λ[p], m])
dt = t[2]-t[1]
integrator = init(prob, algo, saveat=abs.(t[Nₜ÷2+1:-1:1]), dt = dt/10)
# Loop over x and solve the ODE for each initial condition
# Should do this using the ensemble interface in the future
for i in eachindex(x)
# Set up the new initial condition
u0 = ab_dn_t0(x[i])
reinit!(integrator, u0)
# Solve
DiffEqBase.solve!(integrator)
# Multiply by phase factor
rf[Nₜ÷2+1:-1:1, i] .= integrator.sol[1, :].*exp.(+im*(x[i])/4*(2*m-1))
sf[Nₜ÷2+1:-1:1, i] .= integrator.sol[2, :].*exp.(-im*(x[i])/4*(2*m-1))
# Mirror x -> x
rf[Nₜ÷2+2:end, i] .= rf[Nₜ÷2:-1:2, i]
sf[Nₜ÷2+2:end, i] .= sf[Nₜ÷2:-1:2, i]
end
# Compute ψ₁ when p = 1
if p == 1
ψₜ[:,:,n+1] .= ψₜ[:,:,1] .+ (2*(λ[n]' - λ[n]).*sf.*conj.(rf)) ./
(abs2.(rf) + abs2.(sf))
end
else
# Recursion when n != 1
# Calculate r_{n_1, 1} and r_{n-1, p} and similarly for s
r1, s1 = calc_rs(n-1, 1, ψₜ, λ, χ, Ω, Nₓ, Nₜ, x, t, algo)
r2, s2 = calc_rs(n-1, p+1, ψₜ, λ, χ, Ω, Nₓ, Nₜ, x, t, algo)
# Apply the DT equations
rf = ((λ[n-1]' - λ[n-1] ) .* conj.(s1) .* r1 .* s2 +
(λ[p+n-1] - λ[n-1] ) .* abs2.(r1) .* r2 +
(λ[p+n-1] - λ[n-1]') .* abs2.(s1) .* r2 ) ./
(abs.(r1).^2 + abs.(s1).^2);
sf = ((λ[n-1]' - λ[n-1] ) .* s1 .*conj(r1) .* r2 +
(λ[p+n-1] - λ[n-1] ) .* abs2.(s1) .* s2 +
(λ[p+n-1] - λ[n-1]') .* abs2.(r1) .* s2 ) ./
(abs.(r1).^2 + abs.(s1).^2);
if p == 1
# Compute ψₙ when p = 1
ψₜ[:,:,n+1] .= ψₜ[:,:,n] + (2*(λ[n]' - λ[n]).*sf.*conj.(rf)) ./
(abs2.(rf) + abs2.(sf))
end
end
return rf, sf
end
# Solution Parameters
λ = [0.9703563759131811im, 0.9278305490267467im, 0.8522446238645338im, 0.7334225464201684im]
m = 0.0625
χ = [0.7030210896262424, 0.618294940689623, 0.5281884841525464, 0.4276187698960945]
Ω = [0.32801962530412965, 0.6560392506082596, 0.9840588759123894, 1.3120785012165201]
# Simulation Parameters
Nₓ = 1000
xᵣ = -10=>10
dx = (xᵣ.second - xᵣ.first)/Nₓ
x = collect(xᵣ.first:dx:xᵣ.second)[1:end-1]
Nₜ = 1024
N = 3
T = 3*28.7323599983721
dt = T/Nₜ
t = dt * collect((-Nₜ/2:Nₜ/2-1))
algo = Tsit5() <----- Change This
ψ = solve!(λ, m, χ, Ω, Nₓ, Nₜ, x, t, algo)
heatmap(t, x, abs.(ψ)')
```