Consider the following simple SIR model I’ve written in the DifferentialEquations.jl framework:

```
# Import
using DifferentialEquations
# Parameter
𝒫 = [0.1,0.2] # Model Parameters
ℬ = [0.8,0.1, 0.0] # Initial condition
𝒯 = (0.0,365.0) # Time
# Model
function Φ!(du,u,p,t)
S,I,R = u
β,γ = p
du[1] = -β*S*I
du[2] = β*S*I - γ*I
du[3] = γ*I
end
# Problem Definition
problem = ODEProblem(Φ!, ℬ, 𝒯, 𝒫)
# Problem Solution
solution = solve(problem);
```

I would like to extend it with age-stratification in a way that’s both scalable (at least up to 20 age groups per compartment) and compact (vector form).

In the following code chunk I try to define an uncoupled age-specific SIR with only 3 age groups:

```
# Parameter
𝒫 = [0.14,0.12,0.2,0.2] # Model Parameters
ℬ = [0.9,0.9,0.9,0.1,0.1,0.1,0.0,0.0,0.0] # Initial condition
𝒯 = (0.0,365.0) # Time
# Model
function Φ!(du,u,p,t)
s1,s2,s3,i1,i2,i3,r1,r2,r3 = u
β1,β2,β3,γ = p
du[1] = -β1*s1*i1
du[2] = -β2*s2*i2
du[3] = -β3*s3*i3
du[4] = β1*s1*i1-γ*i1
du[5] = β2*s2*i2-γ*i2
du[6] = β3*s3*i3-γ*i3
du[7] = γ*i1
du[8] = γ*i2
du[9] = γ*i3
end
# Problem Definition
problem = ODEProblem(Φ!, ℬ, 𝒯, 𝒫)
# Problem Solution
solution = solve(problem);
```

This works fine but it’s neither compact nor scalable. Can you help me vectorize it in a way to make it so?