# Why the solution of a system of ODEs is different when you run the code a second time?
**URL:** https://discourse.julialang.org/t/why-the-solution-of-a-system-of-odes-is-different-when-you-run-the-code-a-second-time/99709
**Category:** New to Julia
**Tags:** question, numerics, differentialequation
**Created:** [June 1, 2023, 7:36am UTC](https://discourse.julialang.org/t/why-the-solution-of-a-system-of-odes-is-different-when-you-run-the-code-a-second-time/99709 "2023-06-01T07:36:44Z")
**Posts on this page:** 4
**Page:** 1
### Author:  [@atmyo4](https://discourse.julialang.org/u/atmyo4)
#### Post date: [June 1, 2023, 7:36am UTC](https://discourse.julialang.org/t/why-the-solution-of-a-system-of-odes-is-different-when-you-run-the-code-a-second-time/99709/1 "2023-06-01T07:36:44Z")
Hello, I am trying to study the effect of a disturbance on a system of ODEs and when I run the code more than once it gives me a slightly different result compared to when I ran it for the first time. I can’t find the cause of this. When I use callbacks to apply the perturbation instead of solving the systems twice the result is the same as the first run. Can someone please help me figure out what’s going on?
Here is my code, hope it’s not too confusing:
> **Code**
>
> ```julia
> using DifferentialEquations
> using ModelingToolkit
> using Plots
>
> # Define ODE systems
> module model
> using ModelingToolkit
> @parameters mY n gZ1 theta gZ2 bI g1 bC gC
> @variables t Z1(t) Z2(t) X1(t) XC(t)
> D = Differential(t)
>
> @named odeFB = ODESystem([D(Z1) ~ mY - (n * Z1 * Z2) - (gZ1 * Z1),
> D(Z2) ~ (theta * XC) - (n * Z1 * Z2) - (gZ2 * Z2),
> D(X1) ~ (bI * Z1) - (g1 * X1),
> D(XC) ~ (bC * X1) - (gC * XC)])
>
> # External function to obtain xc(t)
> xc(t) = solFBR(t)[4]
> @register_symbolic xc(t)
>
> @named odeNF = ODESystem([D(Z1) ~ mY - (n * Z1 * Z2) - (gZ1 * Z1),
> D(Z2) ~ (theta * xc(t)) - (n * Z1 * Z2) - (gZ2 * Z2),
> D(X1) ~ (bI * Z1) - (g1 * X1),
> D(XC) ~ (bC * X1) - (gC * XC)])
> end
>
> # Set up the parameter values
> @parameters mY n gZ1 theta gZ2 bI g1 bC gC
> p = Dict([
> bI => 0.06, # min-1
> mY => 1 * 60, # min-1 nM
> n => 0.01, # nM-1 min-1
> gZ1 => 0, # nM min-1 [0.0301,0.0408]
> gZ2 => 0, # nM min-1 [0.0301,0.0408]
> theta => 0.3, # min-1
> g1 => 0.1, # min-1
> bC => 1.06*1.05, # min-1 [0.01,0.05565,0.1,0.106]
> gC => 0.1 # nM
> ]);
> pV = collect(pairs(p))
>
> # Set the initial state
> x0 = solve(ODEProblem(model.odeFB,ones(4),1e6,pV); reltol = 1e-12, save_everystep = false).u[2];
>
> # Solve odeFB
> solFB1 = solve(ODEProblem(model.odeFB, x0, 100, pV), alg_hints=[:stiff]; reltol = 1e-12);
> # Pass the solution to the module that contains the differential equations
> @eval model solFBR = Main.solFB1;
> # Solve odeNF
> solNF1 = solve(ODEProblem(model.odeNF, x0, 100, pV), alg_hints=[:stiff]; reltol = 1e-12);
> # Solve odeFB with different initial conditions and tspan
> solFB2 = solve(ODEProblem(model.odeFB, last(solFB1.u), 900, pV), alg_hints=[:stiff]; reltol = 1e-12);
> @eval model solFBR = Main.solFB2;
> # Apply a perturbation to parameter mY
> p[mY] *= 1.05
> pV = collect(pairs(p))
> # Again solve both systems of differential equations
> solFB3 = solve(ODEProblem(model.odeFB, last(solFB1.u), 900, pV), alg_hints=[:stiff]; reltol = 1e-12);
> solNF2 = solve(ODEProblem(model.odeNF, last(solNF1.u), 900, pV), alg_hints=[:stiff]; reltol = 1e-12);
> p[mY] /= 1.05
>
> # Plot
> solFBt = vcat(solFB1.t, solFB3.t .+ 100);
> solFBu = vcat(solFB1[4,:], solFB3[4,:]);
>
> solNFt = vcat(solNF1.t, solNF2.t .+ 100);
> solNFu = vcat(solNF1[4,:], solNF2[4,:]);
>
> plot(solNFt, solNFu, label = "NF")
> plot!(solFBt, solFBu, label = "FB")
>
> ```
Here is the output of the first run:
Here is the output of the second run:
---
### Author:  [@albheim](https://discourse.julialang.org/u/albheim)
#### Post date: [June 1, 2023, 9:22am UTC](https://discourse.julialang.org/t/why-the-solution-of-a-system-of-odes-is-different-when-you-run-the-code-a-second-time/99709/2 "2023-06-01T09:22:57Z")
I tried removing the model module, and then I get consistent results. Not sure why it would have that effect, probably something funky going on with replacing the module and using eval that is over my head…
> **Code**
>
> ```julia
> using DifferentialEquations
> using ModelingToolkit
> using Plots
> using ModelingToolkit
>
> @parameters mY n gZ1 theta gZ2 bI g1 bC gC
> @variables t Z1(t) Z2(t) X1(t) XC(t)
> D = Differential(t)
>
> @named odeFB = ODESystem([D(Z1) ~ mY - (n * Z1 * Z2) - (gZ1 * Z1),
> D(Z2) ~ (theta * XC) - (n * Z1 * Z2) - (gZ2 * Z2),
> D(X1) ~ (bI * Z1) - (g1 * X1),
> D(XC) ~ (bC * X1) - (gC * XC)])
>
> # External function to obtain xc(t)
> xc(t) = solFBR(t)[4]
> @register_symbolic xc(t)
>
> @named odeNF = ODESystem([D(Z1) ~ mY - (n * Z1 * Z2) - (gZ1 * Z1),
> D(Z2) ~ (theta * xc(t)) - (n * Z1 * Z2) - (gZ2 * Z2),
> D(X1) ~ (bI * Z1) - (g1 * X1),
> D(XC) ~ (bC * X1) - (gC * XC)])
>
> # Set up the parameter values
> @parameters mY n gZ1 theta gZ2 bI g1 bC gC
> p = Dict([
> bI => 0.06, # min-1
> mY => 1 * 60, # min-1 nM
> n => 0.01, # nM-1 min-1
> gZ1 => 0, # nM min-1 [0.0301,0.0408]
> gZ2 => 0, # nM min-1 [0.0301,0.0408]
> theta => 0.3, # min-1
> g1 => 0.1, # min-1
> bC => 1.06*1.05, # min-1 [0.01,0.05565,0.1,0.106]
> gC => 0.1 # nM
> ]);
> pV = collect(pairs(p))
>
> # Set the initial state
> x0 = solve(ODEProblem(odeFB,ones(4),1e6,pV); reltol = 1e-12, save_everystep = false).u[2];
>
> # Solve odeFB
> solFB1 = solve(ODEProblem(odeFB, x0, 100, pV), alg_hints=[:stiff]; reltol = 1e-12);
> # Pass the solution to the module that contains the differential equations
> solFBR = solFB1;
> # Solve odeNF
> solNF1 = solve(ODEProblem(odeNF, x0, 100, pV), alg_hints=[:stiff]; reltol = 1e-12);
> # Solve odeFB with different initial conditions and tspan
> solFB2 = solve(ODEProblem(odeFB, last(solFB1.u), 900, pV), alg_hints=[:stiff]; reltol = 1e-12);
> solFBR = solFB2;
> # Apply a perturbation to parameter mY
> p[mY] *= 1.05
> pV = collect(pairs(p))
> # Again solve both systems of differential equations
> solFB3 = solve(ODEProblem(odeFB, last(solFB1.u), 900, pV), alg_hints=[:stiff]; reltol = 1e-12);
> solNF2 = solve(ODEProblem(odeNF, last(solNF1.u), 900, pV), alg_hints=[:stiff]; reltol = 1e-12);
> p[mY] /= 1.05
>
> # Plot
> solFBt = vcat(solFB1.t, solFB3.t .+ 100);
> solFBu = vcat(solFB1[4,:], solFB3[4,:]);
>
> solNFt = vcat(solNF1.t, solNF2.t .+ 100);
> solNFu = vcat(solNF1[4,:], solNF2[4,:]);
>
> plot(solNFt, solNFu, label = "NF")
> plot!(solFBt, solFBu, label = "FB")
>
> ```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [June 1, 2023, 12:00pm UTC](https://discourse.julialang.org/t/why-the-solution-of-a-system-of-odes-is-different-when-you-run-the-code-a-second-time/99709/3 "2023-06-01T12:00:42Z")
> [@atmyo4](#):
>
> ```julia
> # Pass the solution to the module that contains the differential equations
> @eval model solFBR = Main.solFB1;
>
> ```
Yeah this is really funky. Generally if you’re doing something with `@eval` you probably want to consider alternative styles to avoid it.
> [@atmyo4](#):
>
> ```julia
> # Apply a perturbation to parameter mY
> p[mY] *= 1.05
>
> ```
If `p` is the same vector in both problems, then this will mutate both.
---
### Author:  [@atmyo4](https://discourse.julialang.org/u/atmyo4)
#### Post date: [June 2, 2023, 6:13am UTC](https://discourse.julialang.org/t/why-the-solution-of-a-system-of-odes-is-different-when-you-run-the-code-a-second-time/99709/4 "2023-06-02T06:13:58Z")
I think I’ve figured out what’s happening. The first simulation gets interrupted, and since I use the last values of the solution as the initial state for the second simulation, it gives me a different result. In fact, the code did give me a warning, but I ignored it 😛 I suppose I can fix the problem by specifying the algorithm.
> 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. [ODE Solvers · DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Stiff-Problems)).