I’m simulating a model with 8 states using `DifferentialEquations`

and get the expected result:

[Red line is simulation; scatter plot is experimental values.]

Next, I double the model by simulating two identical models (16 states) which are individually independent. When I re-do the simulation, there is a deviation:

These results are from Jupyter notebook. I get the same results if I run the code in VC. I have tried to reduce `reltol`

and `abstol`

, without any effect.

Here is the ODE for the first plot:

```
function mod1!(dx,x,p,t)
S,E,I,C,A,R,CC,xm = x
k0_i,k_e,k_ca,eta,k_r,tm,Tm,N,u = p
k_i = k0_i*xm
k_c = eta*k_ca
k_a = (1-eta)*k_ca
#
dS = -k_i*(I+A)*S/N
dE = k_i*(I+A)*S/N - k_e*E
dI = k_e*E - k_ca*I
dC = k_c*I - k_r*C
dA = k_a*I - k_r*A
dR = k_r*(C+A)
dCC = k_c*I
dxm = (u(t+tm) - xm)/Tm
dx .= [dS,dE,dI,dC,dA,dR,dCC,dxm]
end
```

In this model, `u`

is an `ConstantInterpolation`

function. I call the function as follows:

```
x0 = [S0,E0,I0,C0,A0,R0,CC0,xm0]
#
p = [k0_i,k_e,k_ca,eta,k_r,tm,Tm,N,u]
#
prob = ODEProblem(mod1!,x0,tspan,p)
sol = solve(prob)
```

The model for the second plot is:

```
function mod2!(dx,x,p,t)
S,E,I,C,A,R,CC,xm,S2,E2,I2,C2,A2,R2,CC2,xm2 = x
k0_i,k0_i2,k_e,k_ca,eta,k_r,tm,Tm,tm2,Tm2,N,u,N2,u2 = p
k_i = k0_i*xm
k_i2 = k0_i2*xm2
k_c = eta*k_ca
k_a = (1-eta)*k_ca
#
dS = -k_i*(I+A)*S/N
dE = k_i*(I+A)*S/N - k_e*E
dI = k_e*E - k_ca*I
dC = k_c*I - k_r*C
dA = k_a*I - k_r*A
dR = k_r*(C+A)
dCC = k_c*I
dxm = (u(t+tm) - xm)/Tm
#
dS2 = -k_i2*(I2+A2)*S2/N2
dE2 = k_i2*(I2+A2)*S2/N2 - k_e*E2
dI2 = k_e*E2 - k_ca*I2
dC2 = k_c*I2 - k_r*C2
dA2 = k_a*I2 - k_r*A2
dR2 = k_r*(C2+A2)
dCC2 = k_c*I2
dxm2 = (u2(t+tm2) - xm2)/Tm2
#
dx .= [dS,dE,dI,dC,dA,dR,dCC,dxm,dS2,dE2,dI2,dC2,dA2,dR2,dCC2,dxm2]
end
```

I call the function as follows:

```
x0 = [S0,E0,I0,C0,A0,R0,CC0,xm0,S0,E0,I0,C0,A0,R0,CC0,xm0]
#
p = [k0_i,1.5*k0_i,k_e,k_ca,eta,k_r,tm,Tm,tm,Tm,N,u,N,u]
#
tin = 0.
tfin = 245
tspan = (tin,tfin)
prob = ODEProblem(mod2!,x0,tspan,p)
sol = solve(prob)
```

Observe that in `mod2!`

, the two models (`S`

… `xm`

, and `S2`

… `xm2`

) are (i) identical, and (ii) independent of each other. The only difference when simulating them is that for the second part of `mod2!`

, I change parameter `k0_i`

by multiplying it with `1.5`

.

If I set `k0_i2 = k0_i`

, then I get the correct result for the simulation, but since the models are independent, I should get the same result no matter what `k0_i2`

is??

**Question**: What am I doing wrong??

[I can provide a full set of code, if necessary, but it is rather lengthy.]