# How do I solve a large system of differential equations with indexed functions?

I strongly suspect that matrix methods are the right answer here, but I am having trouble visualizing how to set it up. Imagine a large number of coupled functions:

AB[m,n][t]


where m and n are whole number indices that can range up to hundreds or low thousands, and t is a continuous time variable. Think of AB[m,n] as the concentration of AB of type [m,n], and these concentrations of the different types can evolve with time. In addition, we have a single extra function:

B[t]


that is also a concentration that evolves with time, but for which there is only a single type. Initial conditions are:

AB[a,0][0] = A0
AB[anythingelse][0] = 0
B[0] = B0


Here, a is a constant whole number in the hundreds to low thousands.

AB[m,n] can transform into other types by two processes:

AB[m,n] + B --> AB[m-1,n+b]
AB[m,n] --> AB[m-1,n-1]


Here, b is a constant whole number significantly less than a. This means that the differential equations governing the evolution of the system are:

D[AB[m, n][t], t] == -k1[m, n] AB[m, n][t] B[t]
+k1[m + 1, n - b] AB[m + 1, n - b][t] B[t]
-k2[m, n] AB[m, n][t]
+k2[m + 1, n + 1] AB[m + 1, n + 1][t]
D[B[t], t] == -Sum[k1[m, n] AB[m, n][t],{m,0,a},{n,0,Infinity}] B[t]


Good guesses for the forms for k1 and k2, which we can use for testing, are:

k1[m_] = (m/a) k3
k2[m_, n_] = k4 n/(n + k5/m) + n k6 k1[m]


Here, k3, k4, k5, and k6 are positive real numbers.

How on Earth do I organize this system in a way to numerically solve this system using NDSolve (or ParametricNDSolve) or similar? The ultimate goal will be to have a measurable function that would be something like:

Sum[m AB[m,n],{m,0,a},{n,0,infinity}]


This function would then be fit to experimental data, varying k3-k6, and maybe a and b. Later generalizations may be to have a broader distribution of starting concentrations rather than all being identical, a range of different possible bs, etc.

• Are you really trying to solve a system of ODEs with an infinite number of variables (that is, does n really range from 0 to Infinity)? If not, have you tried creating a matrix equation to represent the ODEs, something like: AB'[t] == m1 . AB[t] . m2 B[t] + n1 . AB[t] . n2, where m1, m2, n1 and n2 encode the information about your k variables. Commented Nov 11, 2019 at 22:43
• Given the mechanism, the maximum value that n could take is a times b. I will look at your suggestion. Commented Nov 11, 2019 at 22:55
• Looks like I found a way to make your suggestion work! Thank you! I will post a completed answer shortly, once I debug one aspect. Commented Nov 12, 2019 at 0:58

a = 20;
b = 4;


Also, my equations above are slightly off. The first transformation should have been:

AB[m,n] + B --> AB[m-1,n+(b-1)]


which changes the differential equations somewhat. First, let's put in our definitions for the rate constants.

k1[m_, n_] = (m/a) k3
k2[m_, n_] = k4 n m/(n m + k5) + n k6 k1[m]


And some initial guess values:

k3 = k4 = k5 = k6 = 1;
AB0 = 10;
B0 = 50;


Now we'll make a vector for all of the species:

ABvector = Flatten@Table[AB[m, n][t], {m, 0, a}, {n, 0, a (b - 1)}];


We can index AB[m,n] via:

index[m_, n_] = m (a (b - 1) + 1) + n + 1;
AB[[index[m,n]]]


Now the k1 and k2 matrices:

k1matrix = SparseArray[
Flatten@Table[{{index[m, n], index[m, n]} -> -k1[m, n],
{index[m - 1, n + (b - 1)], index[m, n]} -> k1[m, n]},
{m, 1, a},
{n, 0, a (b - 1) - (b - 1)}],
{(a + 1) (1 + a (b - 1)), (a + 1) (1 + a (b - 1))}];
k2matrix = SparseArray[
Flatten@Table[{{index[m, n], index[m, n]} -> -k2[m, n],
{index[m - 1, n - 1], index[m, n]} -> k1[m, n]},
{m, 1, a},
{n, 1, a (b - 1)}],
{(a + 1) (1 + a (b - 1)), (a + 1) (1 + a (b - 1))}];


We can now define all of our rate equations (except for B) like this:

rateEq = Thread[D[ABvector, t] == k1matrix.ABvector B[t] + k2matrix.ABvector];


And our initial conditions are (except for B):

initCond = Thread[(ABvector /. t -> 0) == SparseArray[{index[a, 0] -> AB0},
{(a + 1) (1 + a (b - 1))}]];


We can now run NDSolve:

result = NDSolve[Join[rateEq, initCond,
{B[0] == B0, B'[t] == -Sum[k1[m, n] B[t] AB[m, n][t],
{m, 1, a},
{n, 0, a (b - 1) - (b - 1)}]}],
Append[ABvector, B[t]], {t, 0, 10}, WorkingPrecision -> 30];


And we can visualize the results like so:

fn[t_] = Table[AB[m, n][t], {m, 0, a}, {n, 0, a (b - 1)}] /. (result[[1]]);
frameTable = Table[ListPointPlot3D[fn[t], ImageSize -> Large, Filling -> 0,
FillingStyle -> Thickness[0.01], ColorFunction -> "Rainbow",
AxesLabel -> {"n", "m", "Concentration"},
ViewPoint -> {1.3, -2.4, 1.}], {t, 0, 3, 0.1}];
ListAnimate[frameTable]