I am trying to implement a model that is implemented in Matlab in Julia using the DifferentialEquations package. In Julia, the solver stops after a couple of iterations with a warning that an instability has been detected. Here is my (simplified) code:

```
function h(dz, z, p, t)
H, HC, xn = p
# Get value of xn at time t
idx = max(Int(ceil(t/0.05)), 1)
z[:,2:4] = exp.(z[:,2:4])
fv = z[:,3].^(1 ./ H[:,4])
ff = (1 .- (1 .- H[:,5])) .^ (1 ./ z[:,2]) ./ H[:,5]
dz[:,1] = xn[idx] .- H[:,1].*z[:,1] - H[:,2].*(z[:,2] .- 1)
dz[:,2] = z[:,1] ./ z[:,2]
dz[:,3] = (z[:,2] - fv + HC*z[:,5]) ./ (H[:,3].*z[:,3])
dz[:,4] = (ff.*z[:,2] - fv.*z[:,4]./z[:,3] + HC*z[:,6])./(H[:,3].*z[:,4])
dz[:,5] = (-z[:,5] + z[:,3] .- 1)./(H[:,7])
dz[:,6] = (-z[:,6] + z[:,4] .- 1)./(H[:,7])
end
z0 = zeros(2, 6) # initial condition
tspan = (0.0, 600) # time span
H = [0.65 0.41 0.98 0.32 0.34 0.4584 0.5;
0.65 0.41 0.98 0.32 0.34 0.4584 0.5]
HC = [0.0 0.5;
0.0 0.0]
p = [H, HC, xn] # parameters
ts = range(0.05, step=0.05, stop=600) # sample times
ode = ODEProblem(h, z0, tspan, p)
sol = solve(ode, Tsit5(), saveat=ts, tstops=ts)
```

The input parameter xn is the solution of another ODE problem; for simplicity, assume a vector of 0s and 1s:

```
xn = zeros(12000,2)
xn[6000:6004] .= 1
```

The problem occurs with the exponentiation of z[:,2:4] inside the function. (If I remove this line from the function and exponentiate the corresponding elements of the initial states z0, the ODE can be solved. However, that’s not the model I’m trying to implement.) In the first 2 iterations the solutions z are exactly 0.0 and 1.0, then they start fluctuating slightly, until after a few more iterations (about 8) there is a deviation from 1 that is large enough that the values jump on the next iteration and then become infinitely large.

Does anyone have a suggestion what I could do to deal with this instability?