Formulate and solve a system of linear, first-order ODEs in Mathematica

Question

I know this code is wrong, but I am not sure why. I have a system of coupled differential equations. The model is given by the following equations:

$\frac{dS_0} {dt} = (1 - \mu_0) \rho_0 S_0 \frac{S_\gamma}{S_\rho} - \gamma_0 S_0$

$\frac{dS_1} {dt} = \mu_0 \rho_0 S_0 \frac{S_\gamma}{S_\rho} - \gamma_1 S_1 + (1 - \mu_1) \rho_1 S_1 \frac{S_\gamma}{S_\rho}$

$\vdots$

$\frac{dS_n}{dt} = \mu_{n-1} \rho_{n-1} S_{n-1} \frac{S_\gamma}{S_\rho} + (1-\mu_n) \rho_n S_n \frac{S_\gamma}{S_\rho} - \gamma_n S_n$.

Caveat: $\frac{S_\gamma}{S_\rho}$ is a constant term that scales cell birth and cell death rates over all life-stages so the population size remains constant over time. I thought it should be simple to solve the system with Mathematica because I have solved the simplest cases of the ODE system by hand (the 2-stage and 3-stage models) and have an idea that if I put it into matrix form correctly, it will have a triangular structure. Here's my code:

(*Change the recursion limit on Wolfram Mathematica on the Cloud*)
(* Define the number of stages *)
n = 2;
(* Define the state vector *)
S = Array[s, n];
(* Define the coefficient vectors *)
\[Mu] = Array[\[Mu], n];
\[Rho] = Array[\[Rho], n];
\[Gamma] = Array[\[Gamma], n];
(* Define the coefficient matrices *)
A = DiagonalMatrix[-\[Gamma]];
B = DiagonalMatrix[(1 - \[Mu]) \[Rho]];
(* Define the derivative of the state vector *)
dSdt = D[S[t], t];
(* Construct the derivative in terms of matrices A and B *)
ode = dSdt == A.S + B.S \[Gamma]/S \[Rho];
(* Define the initial condition vector *)
S0 = {1, 0}; (* First stage = 1, others = 0 *)
(* Solve the ODE system with the initial condition *)
solution = DSolve[{ode, S[0] == S0}, S[t], t]

If there are n=2 stages, there should be two coupled ODEs; if there are three stages, there should be three coupled ODEs. Is this an inefficient way to input everything to make a matrix ODE system? It seems the way I have formed the matrix ODE system here is either incorrect or missing something because when I run DSolve, it returns it unevaluated. The error message output is: "Function S appears with no arguments." I have searched here a few times to find the best way to formulate similar problems but haven't found anything that jumped out at me as a good way to correct my code.

Answer 1

Starting all over again as correct equations posted.

Since Mathematica index starts from 1 and not zero, it is easier to start your equations from 1 and not 0 as you showed.

Let $\beta=\frac{S_{\gamma}}{S_{\rho}}$ a constant. Then

\begin{align*} S_{1}^{\prime} & =\left( 1-\mu_{1}\right) \rho_{1}S_{1}\beta-\gamma _{1}S_{1}\\ S_{2}^{\prime} & =\mu_{1}\rho_{1}S_{1}\beta-\gamma_{2}S_{2}+\left( 1-\mu _{2}\right) \rho_{2}S_{2}\beta\\ S_{3}^{\prime} & =\mu_{2}\rho_{2}S_{2}\beta-\gamma_{3}S_{3}+\left( 1-\mu _{3}\right) \rho_{3}S_{3}\beta\\ & \vdots\\ S_{n}^{\prime} & =\mu_{n-1}\rho_{n-1}S_{n-1}\beta-\gamma_{n}S_{n}+\left( 1-\mu_{n}\right) \rho_{n}S_{n}\beta \end{align*}

In Matrix form the above becomes

$$ \begin{pmatrix} S_{1}^{\prime}\\ S_{2}^{\prime}\\ S_{3}^{\prime}\\ \vdots\\ S_{n}^{\prime} \end{pmatrix} = \begin{pmatrix} -\gamma_{1}+\left( 1-\mu_{1}\right) \rho_{1}\beta & 0 & 0 & \cdots & \cdots & 0\\ \mu_{1}\rho_{1}\beta & -\gamma_{2}+\left( 1-\mu_{2}\right) \rho_{2}\beta & 0 & \cdots & \cdots & 0\\ 0 & \mu_{2}\rho_{2}\beta & -\gamma_{3}+\left( 1-\mu_{3}\right) \rho_{3}\beta & 0 & \cdots & 0\\ 0 & 0 & \cdots & \cdots & \cdots & 0\\ 0 & 0 & 0 & \cdots & \mu_{n-1}\rho_{n-1}\beta & -\gamma_{n}+\left( 1-\mu _{n}\right) \rho_{n}\beta \end{pmatrix} \begin{pmatrix} S_{1}\\ S_{2}\\ S_{3}\\ \vdots\\ S_{n}% \end{pmatrix} $$

To generate this in Mathematica, a simple loop will do.

n = 2;
S = Array[s, n]
S = Map[#[t] &, S]
dSdt = Map[D[#, t] &, S]
μArray = Array[μ, n]
ρArray = Array[ρ, n]
γArray = Array[γ, n]
eqs = First@Last@Reap@Do[z = (-γArray[[m]] + 
     (1 - μArray[[m]])*ρArray[[m]]*β)*S[[m]]; 
     If[m > 1, z = z + (μArray[[m - 1]]*ρArray[[m - 1]]*β)*S[[m - 1]]]; 
     Sow[dSdt[[m]] == z], {m, 1, n}]

Solve

  DSolve[eqs, S, t]

Here it is for n=3

Solve

 DSolve[eqs, S, t]

This is for n=6

and so on.

old answer

Try this and see if this is what you meant

n = 2;
S = Array[s, n]
S = Map[#[t] &, S]

μArray = Array[μ, n]
ρArray = Array[ρ, n]
γArray = Array[γ, n]

A = DiagonalMatrix[-μArray]
B = DiagonalMatrix[(1 - μArray)  ρArray]
dSdt = Map[D[#, t] &, S]

ode = Thread[dSdt == A . S + B . S  γArray/S  ρArray]

S0 = {1, 0}
ic = MapThread[(#1 == #2) /. t -> 0 &, {S, S0}]

solution = DSolve[{ode, ic}, S, t]