# What is wrong with this system of ODEs?

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:

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

``````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.

``````
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)
``````

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

See this example

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

julia> A
``````

This leads to problems when you acces du like this:

``````du[counter]
``````

``````du[site,stage]
``````

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

``````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)
``````

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