I am trying to solve a second order ODE equation by using SecondOrderODEProblem as problem type. Consistently, I have encountered WARNING: Instability detected. Aborting.

I have tried various solver and still, the problem persists. What are the errors in my code? The same problem is being solved in Matlab without any error.

```
function odefnct(eta,f,df,ddf)
(Phi, iota) = investment(f[2])
# calculate psi
psi_L = eta
psi_R = min(f[2]/df[2] + eta, 1.0)
for n in 1:50
psi = (psi_L + psi_R)/2.0;
amplification = 1.0 - df[2]/f[2]*(psi - eta);
# VOLATILITY COMPUTATION
sigma_eta_eta = sigma*(psi - eta)/amplification; # sigma_eta * eta
sigma_q = sigma_eta_eta*df[2]/f[2];
sigma_theta = sigma_eta_eta*df[1]/f[1];
risk_premium = - sigma_theta*(sigma + sigma_q);
household_premium = (a_ - a)/f[2] + delta - delta_ + risk_premium;
if household_premium > 0.0 # households want to hold more
psi_R = psi;
else
psi_L = psi;
end
end
mu_q = r - (a - iota)/f[2] - Phi + delta - sigma*sigma_q + risk_premium;
mu_eta_eta = -(psi - eta)*(sigma + sigma_q)*(sigma + sigma_q + sigma_theta) + eta*(a - iota)/f[2] + eta*(1.0 - psi)*(delta_ - delta);
qpp = 2.0*(mu_q*f[2] - df[2]*mu_eta_eta)/sigma_eta_eta^2;
thetapp = 2.0*((rho - r)*f[1] - df[1]*mu_eta_eta)/sigma_eta_eta^2;
leverage = psi/eta;
rk = r + risk_premium;
r_k = r + household_premium;
ddf[1] = thetapp
ddf[2] = qpp
end
etaspan = (0.0,1.0)
f0 = [1.0 , q_]
df0 = [-1e+10, 1e+8]
function condition(t,u,integrator)
u[2] - qmax
u[3]
u[4]
end
affect!(integrator) = terminate!(integrator)
cb = ContinuousCallback(condition,affect!)
QL = 0.0 ; QR = 1e+15
global sol
for iter in 1:50
df0[2] = ( QL + QR)/ 2; # this is q'(0)
prob = SecondOrderODEProblem(odefnct, f0, df0, etaspan)
sol = solve(prob, callback = cb)
u1,u2,u3,u4 = sol.u[end]
if abs(u4) < 1e-8 # if q'(eta) has reached zero, we
QL = df0[2]; # increase q'(0)
else # if q(eta) reached qmax or theta'(0) reached 0
QR = df0[2]; # we reduce q'(0)
end
end
```