I’m trying to program a simple matrix-Ricatti ODE solver via the **DifferentialEquations.jl** package and the **RK4** algorithm (I tried the default algorithm as well, to no avail). To benchmark my solver, I’m comparing it to the (user-supplied) MATLAB toolbox “CFH Toolbox”, specifically the function cfaffine.m. The benchmark uses a standard explicit Runge-Kutta 4th-order scheme in a straightforward finite-difference implementation (i.e., no other packages are used). Thus, this would seem like an *almost* apples-to-apples comparison.

However, in my Julia code I start getting errors for initial values with larger imaginary components, something the MATLAB code with the same algorithm handles without issues. Specifically, the errors are “**Warning: dt <= dtmin. Aborting.** There is either an error in your model specification or the true solution is unstable.” These occur regardless of imposing a fixed stepsize (however small or large) or not. Further, even for the intial values that compute without error, as soon as the initial values become complex there is a divergence from the MATLAB benchmark.

What am I misunderstanding about the RK4() implementation? Given that the ODE is numerically straightforward, and especially with a forced fixed stepsize, I don’t understand the divergence from the comparable MATLAB solver or it aborting. Here is my Julia code:

```
# Julia Code
# Ricatti ODE with Runge-Kutta 4th order
# "dt <= dtmin" error occurs with or without fixed stepsize
using DifferentialEquations
using LinearAlgebra
# Original state-space via Duffie-Pan-Singleton 2000 ECMA
K0 = [.02;.01]
K1 = [-1.0 0.0; 0.0 -1.0]
H0 = [0.0 0.0; 0.0 0.0]
H1 = reshape([.16 0.0 0.0 0.0 0.0 0.0 0.0 .09],(2,2,2)) # CUBE
rho0= .005
rho1= [1.0;1.0]
# Quadratic part in Ricatti ODE below
cube_multiply(cube::Array{Float64,3},u) = dropdims(mapslices(i-> u'*i*u, cube; dims=[1,2]);dims=(1,2))
# Ricatti ODE: Stacked system of alpha and beta
function ODE_alphabeta!(du,u,p,t)
du[1] = rho0 - K0'*u[2:end]-1/2*u[2:end]'*H0*u[2:end]; # alpha (singleton dimension)
du[2:end] .= rho1 - K1'*u[2:end] - 1/2*cube_multiply(H1,u[2:end]); # beta
return du
end
function psi(a,x,T)
tspan = (T,0.0) # solving backwards from T to 0
uT = [0;a] # initial values: [alpha(T)=0;beta(T)=a]
prob = ODEProblem(ODE_alphabeta!,uT,tspan)
sol = solve(prob,RK4();dt=0.005,adaptive = false) # using 4th-order Runge-Kutta with same stepsize as MATLAB implementation;
# edited to reflect ChrisRackauckas comment below on adaptive timestep
alpha = sol.u[end][1]
beta = sol.u[end][2:end]
return exp(alpha+beta'*x)
end
# Numerical Example with initial conditions
# Eventual goal: integral \int_0^\infty psi(a0+b0*v*im,x0,T)dv so need solvability for large v
a0 = -[1.0,1.0]
b0 = [1.0,1.0]
x0 = [.02,.01]
T0 = 10
psi(a0,x0,T0) # v=0
psi(a0+5.0*b0*im,x0,T0) # v=5
psi(a0+10.0*b0*im,x0,T0) # v=10
psi(a0+50.0*b0*im,x0,T0) # v=50: this throws a Warning: due to "dt<=dtmin. Aborting."
```

For completeness, below the MATLAB code for replication:

```
% Benchmark MATLAB code (using package "CFH Toolbox", function "cfaffine.m")
K0 = [.02,.01]
K1 = [-1 0 0 -1]
H0 = [0 0 0 0]
H1 = [.16 0 0 0 0 0 0 .09]
R0 = 0.005
R1 = [1 1]
a0 = -[1;1]
b0 = [1;1]
x0 = [.02 .01]
T0 = 10
% Note that the MATLAB solver is by default in the imaginary space
cfaffine(-a0*i,x0,T0,K0,K1,H0,H1,R0,R1,[],[],[],1) % same solution as Julia code above
cfaffine(-a0*i+b0*5,x0,T0,K0,K1,H0,H1,R0,R1,[],[],[],1)
cfaffine(-a0*i+b0*10,x0,T0,K0,K1,H0,H1,R0,R1,[],[],[],1)
cfaffine(-a0*i+b0*50,x0,T0,K0,K1,H0,H1,R0,R1,[],[],[],1)
```