# Instability with DifferentialEquations.jl

**URL:** <https://discourse.julialang.org/t/instability-with-differentialequations-jl/7874>\
**Category:** General Usage\
**Created:** [December 20, 2017, 10:59am UTC](https://discourse.julialang.org/t/instability-with-differentialequations-jl/7874 "2017-12-20T10:59:49Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![newbie](https://avatars.discourse-cdn.com/v4/letter/n/dc4da7/32.png) [@newbie](https://discourse.julialang.org/u/newbie)\
**Post date:** [December 20, 2017, 10:59am UTC](https://discourse.julialang.org/t/instability-with-differentialequations-jl/7874/1 "2017-12-20T10:59:49Z")

</div>

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.

```julia
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

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [December 20, 2017, 2:23pm UTC](https://discourse.julialang.org/t/instability-with-differentialequations-jl/7874/2 "2017-12-20T14:23:15Z")

</div>

What’s `q_`? Also, `investment` isn’t defined.

What method are you using to solve it in MATLAB? Given those initial velocities, is it stiff? Since you didn’t declare a solving method or what type of equation it is, it’s defaulting to an explicit non-stiff solver (like `ode45`). If you’re using something like `ode15s` in MATLAB, you should change the algorithm you’re using in Julia to be something similar (or just pass `alg_hints=[:stiff]`)

---

<div class="post-metadata">

**Author:** ![Rozarspady](https://avatars.discourse-cdn.com/v4/letter/r/ecb155/32.png) [@Rozarspady](https://discourse.julialang.org/u/Rozarspady)\
**Post date:** [October 3, 2020, 2:23pm UTC](https://discourse.julialang.org/t/instability-with-differentialequations-jl/7874/3 "2020-10-03T14:23:01Z")

</div>

#some constants  
const Ip1=12.1299/27.2  
const Ip2=0.78  
const I=7.0e13  
const omg1=0.057  
#const omg2=0.0  
#const ε=0.0  
#const CEP=0.0

const T1=2.0_π/omg1  
 #const T2=2.0_π/omg2  
const E0=sqrt(I/(3.51e16))

---

<div class="post-metadata">

**Author:** ![Rozarspady](https://avatars.discourse-cdn.com/v4/letter/r/ecb155/32.png) [@Rozarspady](https://discourse.julialang.org/u/Rozarspady)\
**Post date:** [October 3, 2020, 2:34pm UTC](https://discourse.julialang.org/t/instability-with-differentialequations-jl/7874/4 "2020-10-03T14:34:43Z")

</div>

I am so sorry that i have made a mistake when I use MATLAB to solve the equation. when I use the same initial condition, I found that ode45() ode15s()in Matlab still cant solve it. When I try to find where the problem is, it seems that when the time t arrive a certain value t0, the step next will go wrong. However, the sol before t0 is Ok,i have calculated du in time t0, there is nothing wrong with them, i think, if du don’t have problems, the iteration is stable ,and i just can’t understand why the solver throw a warning here.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 3, 2020, 6:25pm UTC](https://discourse.julialang.org/t/instability-with-differentialequations-jl/7874/5 "2020-10-03T18:25:04Z")

</div>

This doesn’t seem to be related to the original question, so please make a different post. The issue is certainly in the model somewhere: look for asymptotic properties that could lead to the issue.