# Solve ode problem with parameters described as function of the variables

**URL:** <https://discourse.julialang.org/t/solve-ode-problem-with-parameters-described-as-function-of-the-variables/91042>\
**Category:** New to Julia\
**Tags:** question, ode, differentialequation\
**Created:** [November 30, 2022, 11:29am UTC](https://discourse.julialang.org/t/solve-ode-problem-with-parameters-described-as-function-of-the-variables/91042 "2022-11-30T11:29:20Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![simonluzara](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simonluzara/32/205260_2.png) [@simonluzara](https://discourse.julialang.org/u/simonluzara)\
**Post date:** [November 30, 2022, 11:29am UTC](https://discourse.julialang.org/t/solve-ode-problem-with-parameters-described-as-function-of-the-variables/91042/1 "2022-11-30T11:29:20Z")

</div>

Hello all. I have a sinusoidal excitational model of the hydraulic pumps.  
I have an expression which calculates the density, viscosity and bulk modulus as a fluid(P,T) ,P meaning pressure is a variable.

```julia
function fluid(Pre::Float64,Temp::Float64)
    χ₀=0.001; # Air content percentage in the oil at STP
    Pₛₐₜ=4e6; # Saturation pressure of the oil at stp
    Pᵥₐₚ=1e4; # Vapor pressure of the oil at stp (Shell's data)
    αₜ=7.2e-4; # Volumetric thermal expansion coeff. (/K)
    ρₗ₀=872.0; # oil density at stp
    ρₐ₀=1.225; # air density at stp
    νₚₐ=60.3095; # Kinematic viscosity of liquid oil at operating temperature.(cSt) (Solved by params.m file)
    
    bulk_liq=(10^((0.3766*(log10(νₚₐ)^0.3307))-0.2766)
        +(5.851-0.01382*Temp)*(Pre/1e9))*1e9;
    μₐ=(1.458e-6*(Temp+273.15)^1.5)/(Temp+110.4+273.15); #Sutherland Expression for air
    μₗ₀=(νₚₐ*1e-6)*ρₗ₀*exp(-αₜ*(Temp-T_atm::Float64)); #Dynamic viscosity of liquid oil
    if Pre>=Pₛₐₜ
        rho_f=(ρₗ₀+(χ₀/(1-χ₀))*ρₐ₀)*exp((Pre-p_atm::Float64 )/bulk_liq-αₜ*(Temp-T_atm::Float64));
        Beta_f=bulk_liq;
        mu_f=(νₚₐ*1e-6)*(ρₗ₀+(χ₀/(1-χ₀))*ρₐ₀)*exp((Pre-p_atm::Float64 )/bulk_liq-αₜ*(Temp-T_atm::Float64));
    elseif Pre<Pₛₐₜ && Pre>Pᵥₐₚ
        phi=(Pₛₐₜ-Pre)/(Pₛₐₜ-Pᵥₐₚ);
        rho_f=(ρₗ₀+(χ₀/(1-χ₀))*ρₐ₀)/
            ((phi*(χ₀)/(1-χ₀)*((Temp+273.15)/273.15)*(p_atm::Float64 /Pre)^(1/1.42))+exp((p_atm::Float64 -Pre)/bulk_liq+αₜ*(Temp-T_atm::Float64)));
        
        #Beta_f=(1+(χ₀)/(1-χ₀)*((Temp+273.15)/(273.15))*((p_atm::Float64 /Pre)^(1/1.42))*phi)/
    #(bulk_liq^-1 + (χ₀)/(1-χ₀)*((Temp+273.15)/(273.15))*(p_atm::Float64 /Pre)^(1/1.42)*(1/(Pₛₐₜ-p_atm::Float64 ))*((Pₛₐₜ-Pre)/1.42/Pre+1) );
        
        Beta_f=((phi*(χ₀)/(1-χ₀)*((Temp+273.15)/(273.15))*((p_atm::Float64 /Pre)^(1/1.42)))+exp((p_atm::Float64 -Pre)/bulk_liq+αₜ*(Temp-T_atm::Float64)))/
        (((χ₀)/(1-χ₀)*((Temp+273.15)/(T_atm::Float64 +273.15))*(p_atm::Float64 /Pre)^(1/1.42))*((1/(Pₛₐₜ-Pᵥₐₚ))+(phi/(1.42*Pre)))+
        ((exp((p_atm::Float64 -Pre)/bulk_liq+αₜ*(Temp-T_atm::Float64)))/bulk_liq));
        
        mu_f=((phi*(χ₀)/(1-χ₀)*((Temp+273.15)/273.15)*((p_atm::Float64 /Pre)^(1/1.42))*μₐ)+
        (exp((p_atm::Float64 -Pre)/bulk_liq+αₜ*(Temp-T_atm::Float64))*μₗ₀))/((phi*(χ₀)/(1-χ₀)*((Temp+273.15)/273.15)*((p_atm::Float64 /Pre)^(1/1.42)))+
        exp((p_atm::Float64 -Pre)/bulk_liq+αₜ*(Temp-T_atm::Float64)));
    elseif Pre<=Pᵥₐₚ
        rho_f=(ρₗ₀+(χ₀/(1-χ₀))*ρₐ₀)/
            ((χ₀/(1-χ₀)*((Temp+273.15)/273.15)*((p_atm::Float64 /Pᵥₐₚ)^(1/1.42)))+exp((p_atm::Float64 -Pᵥₐₚ)/bulk_liq+αₜ*(Temp-T_atm::Float64))); 
        Beta_f=1.42*Pre;
        mu_f=μₐ;
    end
    return rho_f,mu_f,Beta_f
end

```

Creating my ode function, I did this:

```julia
function damper!(dPdt,P,pars,t)
    dₚ,dᵣ,Aₚ,Aᵣ,Aₜₚ,hₚ,hᵣ,Δₚ,Δᵣ,lₚₛ,δₚₛ,lᵣₛ,δᵣₛ,
    Hₜ,
    dₒco,dₒre,dₒba,lₒco,lₒre,lₒba,Cfmax_co,Cfmax_re,Cfmax_ba,
    rᵥₐco,lᵥₐco,rₑₓₜco,Aₑₓₜco,kᵥₐco,x₀co,xₛₜᵣco,lₛₚₒco,δco,
    rᵥₐre,lᵥₐre,rₑₓₜre,Aₑₓₜre,kᵥₐre,x₀re,xₛₜᵣre,lₛₚₒre,δre,
    rᵥₐba,lᵥₐba,rₑₓₜba,Aₑₓₜba,kᵥₐba,x₀ba,xₛₜᵣba,lₛₚₒba,δba,
    dₚₐ,lₚₐ,Cfmaxₕ,d_f₁,d_f₂,A_bot,Aₜₒₚ,Cfₗₐₙ,kᵢₙ,x₀,
    V₀co,V₀re,Vg₀,Pg₀,Vres₀,
    amp,freq,=pars;
    
    #--Pressure in the 3 Chambers--#
    Pᵣₑ,Pcₒ,Pᵣₑₛ=P; # Pressure in the rebound chamber,compression chamber and reservoir
    
    #--Imposed Sinusodial Excitation-#
    x=amp*sin(2*pi*freq*t); # Displacement from the middle of the piston 
    u=amp*2*pi*freq*cos(2*pi*freq*t); # Velocity of the piston rod
    
    
    ρ₁,μ₁,β₁=fluid(Pᵣₑ,T);
    ρ₂,μ₂,β₂=fluid(Pcₒ,T);
    ρ₃,μ₃,β₃=fluid(Pᵣₑₛ,T);

    Q_12=ori(Pᵣₑ,Pcₒ,dₒco,lₒco,Cfmax_co,ρ₁,μ₁)+
        ori(Pᵣₑ,Pcₒ,dₒre,lₒre,Cfmax_re,ρ₁,μ₁)+
        sum(spring.(Pᵣₑ,Pcₒ,ρ₁,rᵥₐre,lᵥₐre,rₑₓₜre,Aₑₓₜre,kᵥₐre,x₀re,xₛₜᵣre,lₛₚₒre,δre,3))+
        pi*dₚ*δₚₛ^3/12/lₚₛ/μ₁*abs(Pᵣₑ-Pcₒ)*(Pᵣₑ>=Pcₒ); #Leakage of the piston seals

    Q_21=ori(Pcₒ,Pᵣₑ,dₒco,lₒco,Cfmax_co,ρ₂,μ₂)+
        ori(Pcₒ,Pᵣₑ,dₒre,lₒre,Cfmax_re,ρ₂,μ₂)+
        sum(spring.(Pcₒ,Pᵣₑ,ρ₂,rᵥₐco,lᵥₐco,rₑₓₜco,Aₑₓₜco,kᵥₐco,x₀co,xₛₜᵣco,lₛₚₒco,δco,3))+
        pi*dₚ*δₚₛ^3/12/lₚₛ/μ₂*abs(Pcₒ-Pᵣₑ)*(Pcₒ>=Pᵣₑ); #Leakage of the piston seals

    Q_23=ori(Pcₒ,Pᵣₑₛ,dₒba,lₒba,Cfmax_ba,ρ₂,μ₂)+
        sum(spring.(Pcₒ,Pᵣₑₛ,ρ₂,rᵥₐba,lᵥₐba,rₑₓₜba,Aₑₓₜba,kᵥₐba,x₀ba,xₛₜᵣba,lₛₚₒba,δba,3));

    Q_32=ori(Pᵣₑₛ,Pcₒ,dₒba,lₒba,Cfmax_ba,ρ₃,μ₃)+
        check(Pᵣₑₛ,Pcₒ,Pg₀,dₚₐ,lₚₐ,Cfmaxₕ,d_f₁,d_f₂,A_bot,Aₜₒₚ,Cfₗₐₙ,kᵢₙ,x₀);

    dPdt[1]=β₁/(V₀re-Aₜₚ*x)*(Q_21*(ρ₂/ρ₁)-Q_12
        +Aₜₚ*u);

    dPdt[2]=β₂/(V₀co+Aₚ*x)*(Q_12*(ρ₁/ρ₂)-Q_21+Q_32*(ρ₃/ρ₂)-Q_23
        -Aₚ*u);

    dPdt[3]=(Q_23*(ρ₂/ρ₃)-Q_32)/
        (((Vres₀+Vg₀-((Pg₀/Pᵣₑₛ)^(1/1.42))*Vg₀)/β₃)+(Pg₀^(1/1.42)*Vg₀/1.42/Pᵣₑₛ^(1+1/1.42)));
    
    nothing

    #Force=P_ex*A_tp-P_co*A_p+pi*vis_dyna((P_co+P_ex)/2)*(d_p*h_p/Delta_p+0.5*d_r*h_gu/Delta_r)*abs(u);
end

```

While solving the ode with Rodas5(), I get an error with a huge page, few lines of it being:

```julia
MethodError: no method matching fluid(::ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 3}, ::Float64)
Closest candidates are:
  fluid(::Float64, ::Float64) at In[89]:1

Stacktrace:
  [1] damper!(dPdt::Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 3}}, P::Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 3}}, pars::Vector{Any}, t::Float64)
    @ Main .\In[94]:20
  [2] (::ODEFunction{true, SciMLBase.AutoSpecialize, typeof(damper!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing})(::Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 3}}, ::Vararg{Any})

```

I guess I called the fluid function in a matlab way which would give real values of density, viscosity and bulk modulus. I have checked coupled of tutorials with this kind of a problem where a parameter is function of the variable but I didn’t get a grasp.

Your help will be highly appreciated.

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [November 30, 2022, 11:46am UTC](https://discourse.julialang.org/t/solve-ode-problem-with-parameters-described-as-function-of-the-variables/91042/2 "2022-11-30T11:46:59Z")

</div>

Welcome!

In Julia you don’t need to specify argument types in the method signature.

So in this case I would suggest the following signature

```julia
function fluid(Pre, Temp)

```

---

<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:** [November 30, 2022, 1:09pm UTC](https://discourse.julialang.org/t/solve-ode-problem-with-parameters-described-as-function-of-the-variables/91042/3 "2022-11-30T13:09:07Z")

</div>

what’s happening here is the solver by default tries to use auto differentiation of your code which fails because of the Float64s (which in your case don’t do anything). for models that actually can’t be autodiffed, you can pass autodiff=false to the solver, but here is better to just fix your model.

---

<div class="post-metadata">

**Author:** ![simonluzara](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simonluzara/32/205260_2.png) [@simonluzara](https://discourse.julialang.org/u/simonluzara)\
**Post date:** [November 30, 2022, 1:47pm UTC](https://discourse.julialang.org/t/solve-ode-problem-with-parameters-described-as-function-of-the-variables/91042/4 "2022-11-30T13:47:21Z")

</div>

I read couple of posts on how to optimise a function, and one said I should specify the input variables type.  
I tried to remove all of the `::Float64`, but still the same problem.

---

<div class="post-metadata">

**Author:** ![simonluzara](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simonluzara/32/205260_2.png) [@simonluzara](https://discourse.julialang.org/u/simonluzara)\
**Post date:** [November 30, 2022, 1:49pm UTC](https://discourse.julialang.org/t/solve-ode-problem-with-parameters-described-as-function-of-the-variables/91042/5 "2022-11-30T13:49:15Z")

</div>

It seems yes, but still I don’t understand how does that affect a scalar input when the variable had an explicit value before performing the differentiation.  
I tried to put Rodas5(autodiff=false) and it seems to take ages.

---

<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:** [November 30, 2022, 2:30pm UTC](https://discourse.julialang.org/t/solve-ode-problem-with-parameters-described-as-function-of-the-variables/91042/6 "2022-11-30T14:30:38Z")

</div>

The next problem you’re having is that `T_atm` and `p_atm` are global variables. Pass them in to `fluid` and everything will be a lot faster. I’m pretty sure you mis-interpreted the advice about specifying types.

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [November 30, 2022, 2:32pm UTC](https://discourse.julialang.org/t/solve-ode-problem-with-parameters-described-as-function-of-the-variables/91042/7 "2022-11-30T14:32:53Z")

</div>

> [@simonluzara](#):
>
> I tried to remove all of the `::Float64`, but still the same problem.

Since the error is specifically linked to the method signature, I don’t believe that is correct. You might need to restart julia to override the original method definition (it will take precendence if you don’t).

---

<div class="post-metadata">

**Author:** ![simonluzara](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simonluzara/32/205260_2.png) [@simonluzara](https://discourse.julialang.org/u/simonluzara)\
**Post date:** [November 30, 2022, 2:33pm UTC](https://discourse.julialang.org/t/solve-ode-problem-with-parameters-described-as-function-of-the-variables/91042/8 "2022-11-30T14:33:06Z")

</div>

After performing the calculations with p\_atm and T\_atm as local variables in the `fluid(P)`.

I have received this error about Larger maxiters,

```julia
Warning: Interrupted. Larger maxiters is needed. If you are using an integrator for non-stiff ODEs or an automatic switching algorithm (the default), you may want to consider using a method for stiff equations. See the solver pages for more details (e.g. https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Stiff-Problems).
└ @ SciMLBase C:\Users\490132\.julia\packages\SciMLBase\MVhUm\src\integrator_interface.jl:502

```

I tried CVODE\_BDF, Rodas4 and Vern7 got the same answer.

---

<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:** [November 30, 2022, 2:33pm UTC](https://discourse.julialang.org/t/solve-ode-problem-with-parameters-described-as-function-of-the-variables/91042/9 "2022-11-30T14:33:36Z")

</div>

Also, what matlab solver were you using? It’s possible there’s a better Julia solver for your problem than `Rodas5` (`Rodas5P` is almost always better and depending on tolerances and a bunch of other stuff you might want to be using `QNDF`/`FBDF` or `Rosenbrok23`).

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [November 30, 2022, 2:33pm UTC](https://discourse.julialang.org/t/solve-ode-problem-with-parameters-described-as-function-of-the-variables/91042/10 "2022-11-30T14:33:58Z")

</div>

Is there a reason you are using stiff solvers? I would start with `Tsit5()` if you don’t need a stiff solver

---

<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:** [November 30, 2022, 2:34pm UTC](https://discourse.julialang.org/t/solve-ode-problem-with-parameters-described-as-function-of-the-variables/91042/11 "2022-11-30T14:34:40Z")

</div>

the fact that OP is hitting maxiters implies that the problem is likely stiff.

---

<div class="post-metadata">

**Author:** ![simonluzara](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simonluzara/32/205260_2.png) [@simonluzara](https://discourse.julialang.org/u/simonluzara)\
**Post date:** [November 30, 2022, 2:35pm UTC](https://discourse.julialang.org/t/solve-ode-problem-with-parameters-described-as-function-of-the-variables/91042/12 "2022-11-30T14:35:09Z")

</div>

I used ode23s to start with. It didn’t give me great performance.

---

<div class="post-metadata">

**Author:** ![simonluzara](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simonluzara/32/205260_2.png) [@simonluzara](https://discourse.julialang.org/u/simonluzara)\
**Post date:** [November 30, 2022, 2:36pm UTC](https://discourse.julialang.org/t/solve-ode-problem-with-parameters-described-as-function-of-the-variables/91042/13 "2022-11-30T14:36:28Z")

</div>

The problem is as stiff as a hard bread would be early morning. I can’t tell, I tried with `abstol=1e-8` and `reltol=1e-6`.

---

<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:** [November 30, 2022, 2:42pm UTC](https://discourse.julialang.org/t/solve-ode-problem-with-parameters-described-as-function-of-the-variables/91042/14 "2022-11-30T14:42:05Z")

</div>

What did you use in Matlab?

---

<div class="post-metadata">

**Author:** ![simonluzara](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simonluzara/32/205260_2.png) [@simonluzara](https://discourse.julialang.org/u/simonluzara)\
**Post date:** [November 30, 2022, 2:45pm UTC](https://discourse.julialang.org/t/solve-ode-problem-with-parameters-described-as-function-of-the-variables/91042/15 "2022-11-30T14:45:51Z")

</div>

Method to solve the ode? ode23s and RKF45. The time it takes to calculate especially this check valve equation;

```julia
function check(P_res,P_co,Pg_0,d_o,l_open,Cfmax_o,d_film_1,d_film_2,A_bot,A_top,C_lang,k_in,h_0)
    rho,vis=fluid(P_res/2+P_co/2)[1:2];
    if P_res*A_bot-P_co*A_top-Pg_0*(A_top-A_bot)-k_in*h_0>=0
        Cdv_max=0.27; #Annular Flow maximum value.
        f(h::Float64) = P_res*A_bot - P_co*A_top - k_in*(h+h_0)+
        2*(P_res-P_co)*(C_lang/(pi*d_o^2/4)-A_bot/2*(Cfmax_o*pi*d_o^2/4)^-2)*
        ((Cfmax_o*pi*d_o^2/4)^-2+(Cdv_max*tanh(4*h/vis/100*sqrt(2*(P_res-P_co)*rho))*pi*d_film_1*h)^-2)^-1;
        h_sol=find_zero(f,0.0001,Order2());
        Q_valve=sqrt(2*(P_res-P_co)/rho*
            ((Cfmax_o*pi*d_o^2/4)^-2+(Cdv_max*tanh(4*h_sol/vis/100*sqrt(2*(P_res-P_co)*rho))*pi*d_film_2*h_sol)^-2)^-1);
    else
        Q_valve=0;
    end
end

```

I am sick and tired of this modeling problem. 🙆‍♂️

---

<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:** [November 30, 2022, 3:13pm UTC](https://discourse.julialang.org/t/solve-ode-problem-with-parameters-described-as-function-of-the-variables/91042/16 "2022-11-30T15:13:49Z")

</div>

Rosenbrok23 is the Julia equivalent of ode23s, so it’s probably a decent place to start. How do the Julia and matlab solve times compare? It should be pretty easy to get a fairly substantial speedup while keeping the solvers the same from the language change, but that will require a bit of work to improve your Julia code.

---

<div class="post-metadata">

**Author:** ![simonluzara](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simonluzara/32/205260_2.png) [@simonluzara](https://discourse.julialang.org/u/simonluzara)\
**Post date:** [November 30, 2022, 3:20pm UTC](https://discourse.julialang.org/t/solve-ode-problem-with-parameters-described-as-function-of-the-variables/91042/17 "2022-11-30T15:20:30Z")

</div>

Tried Rosenbrok23 too, same problem with maxiters. I think the problem is the functions that am calling in the ode by `ori(),spring(),check()` where it should have its way of declaring them.  
I am no expert, I have searched everywhere on the internet and I just found one example of parameter as a function of time not a variable.

---

<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:** [November 30, 2022, 3:29pm UTC](https://discourse.julialang.org/t/solve-ode-problem-with-parameters-described-as-function-of-the-variables/91042/18 "2022-11-30T15:29:41Z")

</div>

I recommend going through [How to debug Julia simulation codes (ODEs, optimization, etc.!) - YouTube](https://www.youtube.com/watch?v=g-iOOhh2U6o) to see if any of the debugging techniques there help. There’s a pretty high likelyhood you’ve swapped a sign somewhere that makes the model unstable.

---

<div class="post-metadata">

**Author:** ![simonluzara](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simonluzara/32/205260_2.png) [@simonluzara](https://discourse.julialang.org/u/simonluzara)\
**Post date:** [November 30, 2022, 5:30pm UTC](https://discourse.julialang.org/t/solve-ode-problem-with-parameters-described-as-function-of-the-variables/91042/19 "2022-11-30T17:30:28Z")

</div>

Have watched this JuliaCon and used almost everything being said. I think the problem is the discontinuity by the exogenous functions being called in the damper function. And I can’t make them smooth as they should be.

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [November 30, 2022, 5:38pm UTC](https://discourse.julialang.org/t/solve-ode-problem-with-parameters-described-as-function-of-the-variables/91042/20 "2022-11-30T17:38:00Z")

</div>

Can you predict at which points the exogenous functions are not smooth?

[Next page](https://discourse.julialang.org/t/solve-ode-problem-with-parameters-described-as-function-of-the-variables/91042.md?page=2)