# How to solve instability in ODE problem?

**URL:** <https://discourse.julialang.org/t/how-to-solve-instability-in-ode-problem/73158>\
**Category:** New to Julia\
**Tags:** question, ode, differentialequation\
**Created:** [December 15, 2021, 6:40pm UTC](https://discourse.julialang.org/t/how-to-solve-instability-in-ode-problem/73158 "2021-12-15T18:40:10Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![BirteT](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/birtet/32/31816_2.png) [@BirteT](https://discourse.julialang.org/u/BirteT)\
**Post date:** [December 15, 2021, 6:40pm UTC](https://discourse.julialang.org/t/how-to-solve-instability-in-ode-problem/73158/1 "2021-12-15T18:40:10Z")

</div>

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:

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

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

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [December 15, 2021, 6:59pm UTC](https://discourse.julialang.org/t/how-to-solve-instability-in-ode-problem/73158/2 "2021-12-15T18:59:47Z")

</div>

is this a stiff ODE? If so, you need a stiff solver. Try replacing `Tsit5()` with `AutoTsit5(Rosenbrock23())` which will dynamically switch to a stiff solver if it thinks it needs to. For more details, see [ODE Solvers · DifferentialEquations.jl](https://diffeq.sciml.ai/stable/solvers/ode_solve/)

---

<div class="post-metadata">

**Author:** ![cmarcotte](https://avatars.discourse-cdn.com/v4/letter/c/a3d4f5/32.png) [@cmarcotte](https://discourse.julialang.org/u/cmarcotte)\
**Post date:** [December 15, 2021, 8:58pm UTC](https://discourse.julialang.org/t/how-to-solve-instability-in-ode-problem/73158/3 "2021-12-15T20:58:32Z")

</div>

> [@BirteT](#):
>
> `z[:,2:4] = exp.(z[:,2:4])`

I’m pretty sure your ODE definition should not include something which (selectively?) exponentiates the state values. I think, definitionally, your ODE should be of the form x' = f(x,p,t). In fact, I am not surprised that, since every time `h` is evaluated your state is being (selectively) exponentiated, your state quickly diverges.

In this case, it may help to see the original Matlab code you are re-implementing, or even link to the mathematical model you are using, since I can’t tell from what is included why this exponentiation step is necessary. Hope that is helpful!

---

<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 15, 2021, 9:14pm UTC](https://discourse.julialang.org/t/how-to-solve-instability-in-ode-problem/73158/4 "2021-12-15T21:14:43Z")

</div>

> [@BirteT](#):
>
> `z[:,2:4] = exp.(z[:,2:4])`

Yeah, don’t do that. The only reason why your MATLAB code even works is because this does not do what you expect in MATLAB and it actually ignores this. We don’t ignore your changes, so it will cause bad things. See:

> [@PSA: How to help yourself debug differential equation solving issues](https://discourse.julialang.org/t/psa-how-to-help-yourself-debug-differential-equation-solving-issues/62489):
>
> Debugging differential equation solver issues almost always boils down to doing the same thing, so this is a summary to help you out. For more information on specific issues, check out the FAQ section of the DifferentialEquations.jl documentation which highlights common issues and questions: [https://diffeq.sciml.ai/dev/basics/faq/](https://diffeq.sciml.ai/dev/basics/faq/)How do I debug why the differential equation solver is diverging? dt \<= dtmin. Aborting. There is either an error in your model specification or the true solution i…