Hello, I am calculating a solution for the system using Vern9() in Julia and ode89 in Matlab. I get a message about instability in Julia while there is no such thing in matlab.

If I change the abstol, reltol in julia to 1e-10 and count the trajectory at a time of 250,000 or more, then the trajectory begins to look discrete, this is not observed in matlab.

What could be the reason for this and can it be fixed?

Thanks for the help

Code in Julia:

```
using StaticArrays, DifferentialEquations
function FHN2_try3(u, p ,t)
x1, y1, x2, y2, z= u
ϵ, a, g, k, σ, α, k1, k2 = p
I(ϕ_i) = g * (1.0/(1.0 + exp(k*(cos(σ/2) - cos(ϕ_i - α - σ/2)))))
ρz = k1 + k2 * z ^ 2
ϕ2 = atan(y2, x2)
ϕ1 = atan(y1, x1)
dx1dt = (x1 - x1 ^ 3 / 3 - y1 + I(ϕ2) + ρz * (x2 - x1) ) / ϵ
dy1dt = x1 - a
dx2dt = (x2 - x2 ^ 3 / 3 - y2 + I(ϕ1) + ρz * (x1 - x2) ) / ϵ
dy2dt = x2 - a
dzdt = x1 - x2
return SVector(dx1dt, dy1dt, dx2dt, dy2dt, dzdt)
end
function FHN2_try3_params()
ϵ = 0.01; a = -1.01;
g = 0.1; k = 50.0; σ = 50.0 * pi / 180; α = 160.0 * pi / 180;
k1 = 0.0; k2 = 0.0
return [ ϵ, a, g, k, σ, α, k1, k2]
end
parameters = FHN2_try3_params()
tspan = (0.0, 50_000.0)
parameters[7] = 0.09
parameters[8] = 75.74
u0 = [-0.9859005363852416, -0.635253572091177, -1.0345181027025165, -0.636382088705782, 0.0011285166148596525]
u0 = SVector{5}(u0)
prob = ODEProblem(FHN2_try3, u0, tspan, parameters)
alg = Vern9();
abs_tol = 1e-6;
rel_tol = 1e-3;
max_iters = 1e8;
integrator_setting = (alg = alg, abs_tol = abs_tol, rel_tol = rel_tol, max_iters = max_iters)
sol = solve(prob, integrator_setting.alg, adaptive = true, abstol = integrator_setting.abs_tol, reltol = integrator_setting.rel_tol,
maxiters = integrator_setting.max_iters);
```

Code in Matlab:

```
eps = 0.01; a = -1.01;
g = 0.1; k = 50.0; sigma = 50.0 * pi / 180; alpha = 160.0 * pi / 180;
k1 = 0.09; k2 = 75.74;
tspan = [0 50000];
y0 = [-0.9859005363852416, -0.635253572091177, -1.0345181027025165, -0.636382088705782, 0.0011285166148596525];
opts = odeset('Reltol',1e-3,'AbsTol',1e-6);
[t_solve, solve] = ode89(@(t,yvec) fhn2(t, yvec, eps, a, g, k, sigma, alpha, k1, k2), tspan, y0, opts);
plot(t_solve(length(solve)-10000:length(solve)), solve(length(solve)-10000:length(solve), 1))
xlabel('time')
ylabel('y(t)')
function dx = fhn2(~, v, eps, a, g, k, sigma, alpha, k1, k2)
len = length(v);
x = v(1:2:len - 1);
y = v(2:2:len);
dx = zeros(len, 1);
dx1 = zeros(len, 1);
rho = k1 + k2*v(5)^2;
dx(1:2:len - 1) = (x - x.^3/3 - y + I_vect(atan2(y, x), g, sigma, alpha, k) + rho*[-1, 1; 1, -1]*x)/eps;
dx(2:2:len) = (x - a);
dx(end) = (x(1) - x(2));
end
function res = I_vect(phi, g, sigma, alpha ,k)
G = [0, g;
g, 0];
res = G*(1./(1 + exp(k*(cos(sigma/2) - cos(phi - alpha - sigma/2)))));
end
```