How can I make separate graphs after using NDSolve?

Question

I am solving a system of differential equations with experimental data, I have problem in separating the graphs. Can someone help me?. I attach my code.

It is a system of three differential equations that are connected to each other, however, the parameter "i" causes the total number of differential equations to increase, also, occupies the parameter "i" to read each element of the lists.

The experimental data are grouped in "List".

Thank you.

(*Experimental data*)
p = {0.010121176017058572, 0.010030251430478896, 0.009986915041731468,
    0.009923842643650045, 0.009902780657720081, 0.009939430383633418, 
   0.00998428312848808, 0.010027165278401624, 0.010112291639866421, 
   0.010184341934899909};


(*Experimental data*)
u = {0.252525, 0.249425, 0.247525, 0.24515, 0.243825, 0.243925, 
       0.244225, 0.244475, 0.24575, 0.2467};

(*Experimental data*)
F = {0.0006002300000000001, 0.00067004, 0.00074378, 0.00082196, 
   0.00090455, 0.00099139, 0.0010831500000000002, 0.00117966, 
   0.00128379, 0.00139007};


(*Constants*) 
m = 0.000055;
delta = 33;
gmax = 0.76; 
Oxy = 100;
beta = 12.5;


(*System of differential equations*)
eq = Table[{X[i]'[
      t] == -m p[[i]] F[[i]] Exp[-u[[i]] 0.1] (Y[i][
         t]/(beta + Y[i][t])) (X[i][t] + delta) X[i][t], 
    Y[i]'[t] == - 
        p[[i]] F[[i]] Exp[-u[[i]] 0.1] (Y[i][t]/(beta + Y[i][t])) X[
         i][t] + gmax (1 - (Y[i][t]/Oxy)), 
    Z[i]'[t] == 
     p[[i]] F[[i]] Exp[-u[[i]] 0.1]  (Y[i][t]/(beta + Y[i][t])) X[i][
       t]}, {i, 1, 10}];

(*Initial conditions*)
ini = Table[{X[i][0] == 25, Y[i][0] == 1000, Z[i][0] == 0}, {i, 1, 
    10}];

var = Table[{X[i][t], Y[i][t], Z[i][t]}, {i, 1, 10}];

(*Solution of system of differential equations*)
S = NDSolve[{eq, ini}, var, {t, 0, 1200}];

(*To make the graph*)
Plot[Evaluate[
  Total[Table[{{X[i][t], Y[i][t], Z[i][t]} /. S}, {i, 1, 10}]]], {t, 
  0, 600}]

Answer 1

I tried with DSolve but the RowReduce error tells us that the matrix is not of a regular shape, meaning that the parameters are big and small at the same time.

I got these equations running using NDSolve. To use it you have to set initial conditions which have been added at the end of eq (all equal to zero):

 eq = Flatten[
      Table[{X[i]'[t] == -m p[[i]] F[[i]] Exp[-u[[i]] 0.1] (Y[i][
      t]/(beta + Y[i][t])) (X[i][t] + delta) X[i][t], 
 Y[i]'[t] == -p[[i]] F[[i]] Exp[-u[[i]] 0.1] (Y[i][
       t]/(beta + Y[i][t])) X[i][t] + gmax (1 - (Y[i][t]/Oxy)), 
 Z[i]'[t] == 
  p[[i]] F[[i]] Exp[-u[[i]] 0.1] (Y[i][t]/(beta + Y[i][t])) X[i][
    t], X[i][0] == 0, Y[i][0] == 0, Z[i][0] == 0}, {i, 1, 10}]]

After that I define the variables without the parameter [t], which makes the functions easier to call:

variables = Flatten[Table[{X[i], Y[i], Z[i]}, {i, 1, 10}]]

and then I send it to NDSolve. To solve numerically, you have to set the parameter t over an interval, here from 1 to 10:

solutions = NDSolve[eq, variables, {t, 1, 10}]

You can call each function at a given point (here Y[1] at t=5):

Y[1][5] /. solutions

or plot them,

Plot[Y[1][t] /. solutions, {t, 1, 10}]