Differential Equations: resolving instability in ODE problem

I am trying to transfer a Matlab model to Julia, but am running into instability issues. I have found this discussion:
[Handling Instability When Solving ODE Problems]
and tried some of the suggested actions, but am struggling to diagnose my specific problem. I think that I just lack understanding – I am new to Julia and to working with differential equations – and hope that someone can point me in the right direction.

I have the following (simplified slightly):

using DifferentialEquations, LinearAlgebra, Random, Plots, BenchmarkTools

function f_h!(dz, z, p, t)
    
    xn, PfH, PfHC = p
    
    # Get xn at time t
    dt = 0.05
    idx = max(Int(ceil(t/dt)), 1)
    
    # Exponentiate state variables
    z[:,2:4] = exp.(z[:,2:4])

    # Compute fv and ff
    fv = z[:,3].^(1 ./ H[:,3])
    ff = (1 .- (1 .- H[:,5])) .^ (1 ./ z[:,2]) ./ H[:,5]
    
    # Compute state equations
    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

PfH = [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]
PfHC =  [0.0  -0.693147; 0.0   0.0]

z0 = 0.001*ones(2, 6) # initial condition
tspan = (0.0, 600.0) # time span
p = [xn, PfH, PfHC] # parameters
ts = range(0.05, step=0.05, stop=600.0) # sample times

ode_h = ODEProblem(f_h!, z0, tspan, p)
sol = solve(ode_h, Tsit5(), saveat=ts)
plot(sol)
   

which fails on the first time step with:

┌ Warning: Instability detected. Aborting
└ @ SciMLBase /home/birtet/.julia/packages/SciMLBase/n3U0M/src/integrator_interface.jl:351

The input argument xn is the solution of another ODE Problem (sampled at 0.05 s time steps), which does not suffer from instability issues:

using DifferentialEquations, LinearAlgebra, Random, Plots, BenchmarkTools

function f_n!(x, p, t)
    
    A, B, C, u = p
    
    # Get input at time t
    dt = 0.05
    idx = max(Int(ceil(t/dt)), 1)
    ut = u[idx,:]
    
    # Modulations by inputs
    J = A
    for i in 1:size(B, 3)
        J += B[:,:,i]*ut[i]
    end
    
    J*x + C*ut
end

A = [-1.1  0.0; 0.0 -1.1]
B = zeros(2, 2, 2)
C = [0.3 0.0; 0.0 0.3]

x0 = 0.001*ones(2, 1) # initial condition
tspan = (0.0, 600.0) # time span
p = [A, B, C, u] # parameters
ts = range(0.05, step=0.05, stop=600.0) # sample times

ode_n = ODEProblem(f_n!, x0, tspan, p)
sol = solve(ode_n, Tsit5(), saveat=ts)
xn = sol.u
plot(sol);

(I know that in this particular example the B-matrix seems unnecessary, but its elements aren’t always set to 0.) In f_n! the parameter u represents exogenous/experimental inputs. For reasons of brevity I have not included the code to generate u, but it is a 12000×2 Matrix{Int64} of 0s and 1s.

I have confirmed that the derivative functions in f_h! are exactly the same in MATLAB and Julia. (As a side note, in the MATLAB version f_n! and f_h! are combined into a single ODE problem, but I don’t think that should make a difference - or at least I got the same instability warning when I tried this in Julia.)

When I tried solving ode_h with CVODE_BF(), I got the following error:

[CVODE ERROR]  CVode
  At t = 0 and h = 4.95839e-12, the corrector convergence test failed repeatedly or with |h| = hmin.

When using radau(), there is no error or warning message, but the output is filled with NaNs. Based on the answers to the other post, I think this means that there is an issue with one (or more) of my derivative functions, but I don’t understand what it is. The denominators are all non-zero (at least initially). I’d appreciate any insights into what could be causing the problem.

(I also wanted to read up on defining a Jacobian, but this link no longer works:

http://docs.juliadiffeq.org/latest/features/performance_overloads.html#Declaring-Explicit-Jacobians-1

Ok I found this post:
[differential equations - Julia's DifferentialEquations step size control - Stack Overflow]
and that helped me to figure out where I was going wrong… the solution was to exponentiate the state variables z[:,2:4] outside f_h!