# Nonlinear differential equation struggling to solve

**URL:** <https://discourse.julialang.org/t/nonlinear-differential-equation-struggling-to-solve/123405>\
**Category:** Modelling & Simulations\
**Tags:** differentialequation\
**Created:** [December 3, 2024, 11:12am UTC](https://discourse.julialang.org/t/nonlinear-differential-equation-struggling-to-solve/123405 "2024-12-03T11:12:10Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![Mishal](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mishal/32/214462_2.png) [@Mishal](https://discourse.julialang.org/u/Mishal)\
**Post date:** [December 3, 2024, 11:12am UTC](https://discourse.julialang.org/t/nonlinear-differential-equation-struggling-to-solve/123405/1 "2024-12-03T11:12:10Z")

</div>

Hi Guys

Please refer to the attached code. The code solves a set of coupled differential equations. The code seems to work well for certain input parameters, however sometimes the system seems to be very stiff and does not solve correctly.

I come from a Matlab background, can anyone provide some guidance on how I can alter setting to get better results?

```julia
using DifferentialEquations
using Plots

function statespace!(dx,x,p,t)
    l,m,g,I,kc,kb,C1,C2 = p

    A = x[2]^2*l*sin(x[1])+(-kb*x[3]-C1*x[4])/m
    B = -kc*x[1]-m*g*l*sin(x[1])-C2*x[2]
    C = (1-(m*l^2*cos(x[1])^2)/(l^2*m+I))

    theta_dotdot = ((-m*A*l*cos(x[1]+B)/(l^2*m+I)))*1/C

    x_dotdot = x[2]^2*l*sin(x[1])-theta_dotdot*l*cos(x[1])+(-kb*x[3]-C1*x[4])/m

    dx[1] = x[2]
    dx[2] = theta_dotdot
    dx[3] = x[4]
    dx[4] = x_dotdot  

end

m = 230
l = 5
g = 9.81
I = 230*(l/2)^2
kc = 795774.71
kb = 236*1000
D1 = 0.04*(m+kb)*0
D2 = 0.04*(m+kc)*0

x0 = [0 , 0.0, 0.0, 2.5]  
tspan = (0.0,1) # time span

p = (m,l,g,I,kc,kb,D1,D2)

prob = ODEProblem(statespace!, x0, tspan, p)
sol = solve(prob,Rodas5(),saveat = 0.001,reltol=1e-8,abstol=1e-8)

plt_1 = plot()
plot!(plt_1, sol.t,sol[1,:],xlabel = "Time (sec)",ylabel = "Angular displacement")
display(plot(plt_1))

plt_2 = plot()
plot!(plt_2, sol.t,sol[2,:],xlabel = "Time (sec)",ylabel = "Angular velocity")
display(plot(plt_2))

plt_3 = plot()
plot!(plt_3, sol.t,sol[3,:],xlabel = "Time (sec)",ylabel = "Linear displacement")
display(plot(plt_3))

plt_4 = plot()
plot!(plt_4, sol.t,sol[4,:],xlabel = "Time (sec)",ylabel = "Linear velocity")
display(plot(plt_4))

```

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [December 3, 2024, 1:29pm UTC](https://discourse.julialang.org/t/nonlinear-differential-equation-struggling-to-solve/123405/2 "2024-12-03T13:29:59Z")

</div>

The simplest thing to try is to experiment with solvers, e.g., a high-order one for stiff systems like `RadauIIA9()`:

```julia
sol = solve(prob, RadauIIA9(), reltol=1e-8, abstol=1e-8)

```

You may also try `FBDF, QNDF` and the other Rodas solvers like `Rodas5P()`