I am trying to solve a set of coupled nonlinear equations, and the system has physical instability. The equations of motion are

i\frac{d\rho(\cos\theta,t)}{dt}=[H,\rho(\cos\theta,t)]

i\frac{d\bar{\rho}(\cos\theta,t)}{dt}=[\bar{H},\bar{\rho}(\cos\theta,t)]

where H and \bar{H} are functions of \rho and \bar{\rho} making the system highly nonlinear and the system has a physical instability. The solution grows exponentially. I solve it by discretizing over \cos\theta and I end up with a few thousand odes. If I solve it using KenCarp47 the solution is stable. Are these algorithms mistaking a real instability for a numerical one and avoiding it somehow? I have pasted by Julia code below.

```
using DifferentialEquations
using OrdinaryDiffEq
using Plots
using Sundials
using LinearSolve
const Î¸ = 1.0e-6
const Energy=50
const Î”msq = 2.5e-3*(1.0e6/197.3269631)
const Hvac=(Î”msq/(4.0*Energy))*[[-cos(2.0*Î¸)+0im, sin(2.0*Î¸)+0im] [sin(2.0*Î¸)+0im, cos(2.0*Î¸)+0im]]
const c = 3e5
const mup = 1.0e5
@time function eoms!(du,u,p,t)
println(t)
AngleBins=p[1]
for i = 1:AngleBins*8
du[i]=0.0e0
end
Hself0 = zeros(ComplexF64,2,2)
Hself1=zeros(ComplexF64,2,2)
dcosth=2.0e0/AngleBins
costharr=range(-1.0+dcosth/2.0e0,1.0-dcosth/2.0e0,AngleBins)
for i = 0:AngleBins-1
rho=[[u[i*8+1]+0im, u[i*8+3]-1im*u[i*8+4]] [u[i*8+3]+1im*u[i*8+4],u[i*8+2]+0im]]
rhobar=[[u[i*8+5]+0im, u[i*8+7]-1im*u[i*8+8]] [u[i*8+7]+1im*u[i*8+8],u[i*8+6]+0im]]
Hself0=Hself0+(rho-rhobar)*dcosth
Hself1=Hself1+(rho-rhobar)*costharr[i+1]*dcosth
end
for i = 0:AngleBins-1
rho=[[u[i*8+1]+0im, u[i*8+3]-1im*u[i*8+4]] [u[i*8+3]+1im*u[i*8+4],u[i*8+2]+0im]]
rhobar=[[u[i*8+5]+0im, u[i*8+7]-1im*u[i*8+8]] [u[i*8+7]+1im*u[i*8+8],u[i*8+6]+0im]]
Hself=mup*(Hself0-Hself1*costharr[i+1])
dudt=-c*1im*((Hvac+Hself)*rho-rho*(Hvac+Hself))
dudtbar=-c*1im*((-1.0e0*Hvac+Hself)*rhobar-rhobar*(-1.0e0*Hvac+Hself))
du[i*8+1]=real(dudt[1,1])
du[i*8+2]=real(dudt[2,2])
du[i*8+3]=real(dudt[1,2])
du[i*8+4]=imag(dudt[1,2])
du[i*8+5]=real(dudtbar[1,1])
du[i*8+6]=real(dudtbar[2,2])
du[i*8+7]=real(dudtbar[1,2])
du[i*8+8]=imag(dudtbar[1,2])
end
end
AngleBins=100
size=8*AngleBins
dcosth=2.0e0/AngleBins
costharr=range(-1.0+dcosth/2.0e0,1.0-dcosth/2.0e0,AngleBins)
u0=zeros(Float64,size)
for i=0:AngleBins-1
u0[8*i+1]=0.5e0
u0[8*i+5]=0.47e0+0.05*exp(-(costharr[i+1]-1.0e0)*(costharr[i+1]-1.0e0))
end
p=[AngleBins]
tspan = (0.0,1e-6)
prob = ODEProblem(eoms!,u0,tspan,p)
#reltolp=1e-10
#abstolp=1e-10
#@time sol1 = solve(prob,ARKODE(Sundials.Explicit(),etable = Sundials.DORMAND_PRINCE_7_4_5),reltol=reltolp,abstol=abstolp,saveat=1e-8)
#reltolp=1e-10
#abstolp=1e-10
#@time sol2 = solve(prob,QNDF(),reltol=reltolp,abstol=abstolp,saveat=1e-8)
reltolp=1e-12
abstolp=1e-12
@time sol1 = solve(prob,CVODE_Adams(linear_solver=:GMRES,krylov_dim=30),reltol=reltolp,abstol=abstolp,saveat=1e-8)
@time sol3 = solve(prob,KenCarp47(linsolve=KrylovJL_GMRES(),concrete_jac=true),reltol=reltolp,abstol=abstolp,saveat=1e-8)
#@time sol1 = solve(prob,TRBDF2(linsolve=KrylovJL_GMRES(),concrete_jac=true),reltol=reltolp,abstol=abstolp,saveat=1e-8)
using Plots
plot(sol1,vars=(2),yscale=:log10,ylims=(1e-15,1e-3),label="CVODE_Adams")
plot!(sol3,vars=(2),yscale=:log10,ylims=(1e-15,1e-3),label="KenCarp47")
```