# What is wrong with this system of ODEs?

**URL:** <https://discourse.julialang.org/t/what-is-wrong-with-this-system-of-odes/93924>\
**Category:** Specific Domains\
**Tags:** diffeq, ode\
**Created:** [February 2, 2023, 2:07pm UTC](https://discourse.julialang.org/t/what-is-wrong-with-this-system-of-odes/93924 "2023-02-02T14:07:51Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ali\_Vahdati](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ali_vahdati/32/38715_2.png) [@Ali\_Vahdati](https://discourse.julialang.org/u/Ali_Vahdati)\
**Post date:** [February 2, 2023, 2:07pm UTC](https://discourse.julialang.org/t/what-is-wrong-with-this-system-of-odes/93924/1 "2023-02-02T14:07:51Z")

</div>

I have a system with 5 linear ODEs. The code runs and the results are just what I expect. In a second system, I expand the first system by duplicating the five equations. There are now 10 equations where there is no connection between the first and the second five equations.  
I expected that the the second system gives me identical results to the first system, but duplicated. Instead, I get the following error:

```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://diffeq.sciml.ai/stable/solvers/ode_solve/#Stiff-Problems).

└ @ SciMLBase C:\Users\arvka\.julia\packages\SciMLBase\OfgnG\src\integrator_interface.jl:491

retcode: MaxIters

Interpolation: specialized 4th order "free" interpolation, specialized 2nd order "free" stiffness-aware interpolation

t: 999977-element Vector{Float64}:

```

I tried larger maxiters (1e7) but it didn’t help.

Here is the first system:

```julia
using Pkg
Pkg.activate(".")
using LinearAlgebra
using DifferentialEquations

u0 = [1800.0 1800.0 1800.0 1800.0 1800.0]
p = Any[[5.846581865622962, 0.998611, 0.971057, 0.683772, 0.00629], [0.001, 0.001, 0.001, 0.001, 0.000322], 64000]

function model1(du, u, p, t)
  nsites = 1

  r, d, K = p

  counter = 0
  for site in 1:nsites

    # Stages 1
    counter += 1
    stage = 1
    prev_stage = 5

    du[counter] = (r[stage] * u[site, prev_stage]) -
                  (r[stage+1] * u[site, stage]) -
                  (d[stage] * u[site, stage])

    # Stage 2
    counter += 1
    stage = 2
    prev_stage = 1
    du[counter] = (r[stage] * u[site, prev_stage]) -
                  (r[stage+1] * u[site, stage]) -
                  (d[stage] * u[site, stage])

    #stage 3 
    counter += 1
    stage = 3
    prev_stage = 2
    du[counter] = (r[stage] * u[site, prev_stage] * ((K - u[site, stage]) / K)) -
                  (r[stage+1] * u[site, stage]) -
                  (d[stage] * u[site, stage])

    # stage 4 
    counter += 1
    stage = 4
    prev_stage = 3
    du[counter] = (r[stage] * u[site, prev_stage] * ((K - u[site, stage]) / K)) -
                  (r[stage+1] * u[site, stage]) -
                  (d[stage] * u[site, stage])

    # stage 5 
    counter += 1
    stage = 5
    prev_stage = 4
    du[counter] = (r[stage] * u[site, prev_stage] * ((K - u[site, stage]) / K)) -
                  (0.0 * u[site, stage]) -
                  (d[stage] * u[site, stage])

  end
end

tspan = (1.0, 200.0)
prob = ODEProblem(model1, u0, tspan, p)
sol = solve(prob)

```

And here is the second system. Note that the only difference in this case is `nsites=2`, instead of 1 and u0 is a matrix with two rows.

```julia

u02 = [1800.0 1800.0 1800.0 1800.0 1800.0; 1800.0 1800.0 1800.0 1800.0 1800.0]
p2 = [[5.846581865622962, 0.998611, 0.971057, 0.683772, 0.00629], [0.001, 0.001, 0.001, 0.001, 0.000322], 64000.0]

function model2(du, u, p, t)
  nsites = 2

  r, d, K = p

  counter = 0
  for site in 1:nsites

    # Stages 1
    counter += 1
    stage = 1
    prev_stage = 5

    du[counter] = (r[stage] * u[site, prev_stage]) -
                  (r[stage+1] * u[site, stage]) -
                  (d[stage] * u[site, stage])

    # Stage 2
    counter += 1
    stage = 2
    prev_stage = 1
    du[counter] = (r[stage] * u[site, prev_stage]) -
                  (r[stage+1] * u[site, stage]) -
                  (d[stage] * u[site, stage])

    #stage 3 
    counter += 1
    stage = 3
    prev_stage = 2
    du[counter] = (r[stage] * u[site, prev_stage] * ((K - u[site, stage]) / K)) -
                  (r[stage+1] * u[site, stage]) -
                  (d[stage] * u[site, stage])

    # stage 4 
    counter += 1
    stage = 4
    prev_stage = 3
    du[counter] = (r[stage] * u[site, prev_stage] * ((K - u[site, stage]) / K)) -
                  (r[stage+1] * u[site, stage]) -
                  (d[stage] * u[site, stage])

    # stage 5 
    counter += 1
    stage = 5
    prev_stage = 4
    du[counter] = (r[stage] * u[site, prev_stage] * ((K - u[site, stage]) / K)) -
                  (0.0 * u[site, stage]) -
                  (d[stage] * u[site, stage])

  end
end

tspan = (1.0, 200.0)
prob2 = ODEProblem(model2, u02, tspan, p2)
sol = solve(prob2)

```

---

<div class="post-metadata">

**Author:** ![halejulia](https://avatars.discourse-cdn.com/v4/letter/h/fbc32d/32.png) [@halejulia](https://discourse.julialang.org/u/halejulia)\
**Post date:** [February 2, 2023, 4:33pm UTC](https://discourse.julialang.org/t/what-is-wrong-with-this-system-of-odes/93924/2 "2023-02-02T16:33:55Z")

</div>

I think the problem is how you access du or how you defined u02:

See this example

```julia
julia> A = [1 2 3; 4 5 6]
2×3 Matrix{Int64}:
 1 2 3
 4 5 6

julia> A[2]
4 # instead of 2

```

This leads to problems when you acces du like this:

```julia
du[counter]

```

you could use this instead.

```julia
du[site,stage]

```

Here is the full code with nsites=2 that works correctly I think.

```julia
u02 = [1800.0 1800.0 1800.0 1800.0 1800.0; 1800.0 1800.0 1800.0 1800.0 1800.0]
p2 = [[5.846581865622962, 0.998611, 0.971057, 0.683772, 0.00629], [0.001, 0.001, 0.001, 0.001, 0.000322], 64000.0]

function model2(du, u, p, t)
  nsites = 2

  r, d, K = p

  counter = 0
  for site in 1:nsites

    # Stages 1
    counter += 1
    stage = 1
    prev_stage = 5

    du[site,stage] = (r[stage] * u[site, prev_stage]) -
                  (r[stage+1] * u[site, stage]) -
                  (d[stage] * u[site, stage])

    # Stage 2
    counter += 1
    stage = 2
    prev_stage = 1
    du[site,stage] = (r[stage] * u[site, prev_stage]) -
                  (r[stage+1] * u[site, stage]) -
                  (d[stage] * u[site, stage])

    #stage 3 
    counter += 1
    stage = 3
    prev_stage = 2
    du[site,stage] = (r[stage] * u[site, prev_stage] * ((K - u[site, stage]) / K)) -
                  (r[stage+1] * u[site, stage]) -
                  (d[stage] * u[site, stage])

    # stage 4 
    counter += 1
    stage = 4
    prev_stage = 3
    du[site,stage] = (r[stage] * u[site, prev_stage] * ((K - u[site, stage]) / K)) -
                  (r[stage+1] * u[site, stage]) -
                  (d[stage] * u[site, stage])

    # stage 5 
    counter += 1
    stage = 5
    prev_stage = 4
    du[site,stage] = (r[stage] * u[site, prev_stage] * ((K - u[site, stage]) / K)) -
                  (0.0 * u[site, stage]) -
                  (d[stage] * u[site, stage])

  end
end

tspan = (1.0, 200.0)
prob2 = ODEProblem(model2, u02, tspan, p2)
sol = solve(prob2)

```

---

<div class="post-metadata">

**Author:** ![Ali\_Vahdati](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ali_vahdati/32/38715_2.png) [@Ali\_Vahdati](https://discourse.julialang.org/u/Ali_Vahdati)\
**Post date:** [February 2, 2023, 7:12pm UTC](https://discourse.julialang.org/t/what-is-wrong-with-this-system-of-odes/93924/3 "2023-02-02T19:12:14Z")

</div>

Oh I didn’t know that du has the same shape as u. This solved the problem. Thank you!