How to speed up NDSolve when there are ~400k equations?

I tested the code as it is without optimization. I give a report and general recommendations: Test 1

In[11]:= Equations[tMin, tMax, 200, 20]; // AbsoluteTiming Out[11]= {7.76804, Null} 

Test 2

In[12]:= Equations[tMin, tMax, 200, 200]; // AbsoluteTiming Out[12]= {86.8025, Null} 

Test 3

In[13]:= Equations[tMin, tMax, 2000, 20]; // AbsoluteTiming Out[13]= {109.039, Null} 

Test 4

In[14]:= Equations[tMin, tMax, 2000, 200]; // AbsoluteTiming Out[14]= {2174.46, Null} 

And so it took 36 minutes, 30-50% of the CPU and 12GB of memory (maximum). Most of the time 1 kernel worked. I think that we can reduce the time to 3-5 minutes. Version report:

 $Version Out[]= "12.0.0 for Microsoft Windows (64-bit) (April 6, 2019)" 

Let's start the optimization. Here we are dealing with the Riccati matrix equation $$c'=c.A.c+B1.c+c.B2...(1) $$ where c is a 2000x200 matrix. Therefore we used Method -> {"EquationSimplification" -> "Residual"to solve it. Then it takes about 36 min on my ASUS laptop. We can reduce time by simple replacement t->tmax t. The corresponding code is

PrepareEquations[tMin_, tMax_, NoC_, NoN_] := Module[{i, j,(*C and N species*)\[Alpha]C = 1,(*C monomer*)\[Alpha]N = 1,(*N monomer*)nCmin = 0, nCmax = NoC,(*Number of C atoms*) nNmin = 0, nNmax = NoN,(*Number of N atoms*) c0 = 10^-30(*For numerical stability,avoid null value*)}, Print["Function: PrepareEquations"]; Clear[c, a, b, F, Fc, Fn]; Cinit = 0.00439673;(*Initial C concentration*) Ninit = 0.00761477;(*Initial N concentration*) k = 1.380649*10^-23;(*Boltzmann constant,unit J/K*) T = 723;(*Temperature,unit K*) R[nC_, nN_] := ((3*(nC + nN)*0.0118)/(4 \[Pi]))^(1/3);(*Radius*) D\[Alpha][i_, j_] := 0;(*Diffusion coefficients*) D\[Alpha][\[Alpha]C, 0] := 8.0*10^-4; D\[Alpha][0, \[Alpha]N] := 1.0*10^-4; k\[Alpha][nC_, nN_, i_, j_] := 4 \[Pi]*(R[nC, nN] + R[i, j])/0.0118*D\[Alpha][i, j]; (*Calculation of energy*)RsC = 0.02780299; \[Sigma]C = 4.09324*10^-19; Fc[nC_] := 4 \[Pi]*(R[nC, 0])^2*\[Sigma]C*(1 - RsC/R[nC, 0]); \[Sigma]N = 5.54492*10^-20; RsN = -0.14758424; Fn[nN_] := 4 \[Pi]*(R[0, nN])^2*\[Sigma]N*(1 - RsN/R[0, nN]); F[nC_, nN_] := -k*T* Log[(nC + nN)!/(nC!*nN!)] + (nC*Fc[nC + nN] + nN*Fn[nC + nN])/(nC + nN); b[nC_, nN_, i_, j_] := k\[Alpha][nC, nN, i, j]*c[i, j][t];(*Condensation rate*) a[nC_, nN_, i_, j_] := k\[Alpha][nC - i, nN - j, i, j]* Exp[-((F[nC - i, nN - j] + F[i, j] - F[nC, nN])/(k* T))];(*Emission rate*)a[i_, j_, i_, j_] := 0;(*Monomer*) a[0, nN_, 1, j_] := 0; a[nC_, 0, i_, 1] := 0; Clear[CL, Eqk, EqD, Eq]; With[{},(*Initial conditions*)CL[nC_, nN_] := c[nC, nN][tMin] == 0; (*Definition of the equilibrium equations*) EqD[nC_, nN_] := c[nC, nN]'[t]; Eqk[nC_, nN_] := Sum[b[nC - i, nN - (1 - i), i, 1 - i]* c[nC - i, nN - (1 - i)][ t] - (b[nC, nN, i, 1 - i] + a[nC, nN, i, 1 - i])* c[nC, nN][t] + a[nC + i, nN + (1 - i), i, 1 - i]* c[nC + i, nN + (1 - i)][t], {i, {0, 1}}]; Eq[nC_, nN_] := EqD[nC, nN] == tMax Eqk[nC, nN]; (*Definition of equilibrium concentrations for each special size*) Module[{nC, nN}, nC = nCmax; For[nN = 1, nN < nNmax, nN++, Eqk[nC, nN] = (b[nC, nN - 1, 0, 1]* c[nC, nN - 1][t] - (b[nC, nN, 0, 1] + a[nC, nN, 0, 1])* c[nC, nN][t] + a[nC, nN + 1, 0, 1]* c[nC, nN + 1][t]) + (b[nC - 1, nN, 1, 0]*c[nC - 1, nN][t] - a[nC, nN, 1, 0]*c[nC, nN][t])]; nN = nNmax; For[nC = 1, nC < nCmax, nC++, Eqk[nC, nN] = (b[nC - 1, nN, 1, 0]* c[nC - 1, nN][t] - (b[nC, nN, 1, 0] + a[nC, nN, 1, 0])* c[nC, nN][t] + a[nC + 1, nN, 1, 0]* c[nC + 1, nN][t]) + (b[nC, nN - 1, 0, 1]*c[nC, nN - 1][t] - a[nC, nN, 0, 1]*c[nC, nN][t])]; nC = 0; For[nN = 2, nN < nNmax, nN++, Eqk[nC, nN] = (b[nC, nN - 1, 0, 1]* c[nC, nN - 1][t] - (b[nC, nN, 0, 1] + a[nC, nN, 0, 1])* c[nC, nN][t] + a[nC, nN + 1, 0, 1]* c[nC, nN + 1][t]) + (-b[nC, nN, 1, 0] c[nC, nN][t] + a[nC + 1, nN, 1, 0]*c[nC + 1, nN][t])]; nN = 0; For[nC = 2, nC < nCmax, nC++, Eqk[nC, nN] = (b[nC - 1, nN, 1, 0]* c[nC - 1, nN][t] - (b[nC, nN, 1, 0] + a[nC, nN, 1, 0])* c[nC, nN][t] + a[nC + 1, nN, 1, 0]* c[nC + 1, nN][t]) + (-b[nC, nN, 0, 1] c[nC, nN][t] + a[nC, nN + 1, 0, 1]*c[nC, nN + 1][t])]; nC = nCmax; nN = 0; Eqk[nC, nN] = (-b[nC, nN, 0, 1]*c[nC, nN][t] + a[nC, nN + 1, 0, 1]*c[nC, nN + 1][t]) + (b[nC - 1, nN, 1, 0]* c[nC - 1, nN][t] - a[nC, nN, 1, 0]*c[nC, nN][t]); nC = 0; nN = nNmax; Eqk[nC, nN] = (-b[nC, nN, 1, 0]*c[nC, nN][t] + a[nC + 1, nN, 1, 0]*c[nC + 1, nN][t]) + (b[nC, nN - 1, 0, 1]* c[nC, nN - 1][t] - a[nC, nN, 0, 1]*c[nC, nN][t]); nC = nCmax; nN = nNmax; Eqk[nC, nN] = (b[nC, nN - 1, 0, 1]*c[nC, nN - 1][t] - a[nC, nN, 0, 1]*c[nC, nN][t]) + (b[nC - 1, nN, 1, 0]* c[nC - 1, nN][t] - a[nC, nN, 1, 0]*c[nC, nN][t]); nC = 0; nN = 0; CL[nC, nN] = c[nC, nN][tMin] == 0; Eqk[nC, nN] = 0; nC = \[Alpha]C; nN = 0; CL[nC, nN] = c[nC, nN][tMin] == Cinit + c0; Eqk[nC, nN] = (-b[nC, nN, 1, 0]*c[nC, nN][t] + a[nC + 1, nN, 1, 0]*c[nC + 1, nN][t]) + (-b[nC, nN, 0, 1]* c[nC, nN][t] + a[nC, nN + 1, 0, 1]* c[nC, nN + 1][ t]) - (Sum[(b[x, y, nC, nN] - a[x, y, nC, nN])* c[x, y][t], {x, 0, nCmax - 1}, {y, 0, nNmax}] + Sum[(-a[nCmax, i, nC, nN]*c[nCmax, i][t]), {i, 0, nNmax}]); nC = 0; nN = \[Alpha]N; CL[nC, nN] = c[nC, nN][tMin] == Ninit + c0; Eqk[nC, nN] = (-b[nC, nN, 1, 0]*c[nC, nN][t] + a[nC + 1, nN, 1, 0]*c[nC + 1, nN][t]) + (-b[nC, nN, 0, 1]* c[nC, nN][t] + a[nC, nN + 1, 0, 1]* c[nC, nN + 1][ t]) - (Sum[(b[x, y, nC, nN] - a[x, y, nC, nN])* c[x, y][t], {x, 0, nCmax}, {y, 0, nNmax - 1}] + Sum[(-a[i, nNmax, nC, nN]*c[i, nNmax][t]), {i, 0, nCmax}]);]; Bc := Flatten[ParallelTable[ c[nC, nN], {nC, 0, nCmax}, {nN, 0, nNmax}]] // Evaluate]; DistributeDefinitions[Eq, CL]; Sys = Join[ Flatten[ParallelTable[ Eq[nC, nN], {nC, 0, nCmax}, {nN, 0, nNmax}]], Flatten[ParallelTable[ CL[nC, nN], {nC, 0, nCmax}, {nN, 0, nNmax}]]] // Flatten // Evaluate; {Sys, Bc}]; tMin = 0; tMax = 7.75*10^5; 

Here we have replaced only one equation Eq[nC_, nN_] := EqD[nC, nN] == tMax Eqk[nC, nN];.Now we call

PrepareEquations[tMin, tMax, 2000, 200]; // AbsoluteTiming 

And it takes {536.24, Null}. Now we solve the system on {t, 0, 1}:

RSol = NDSolve[Sys, Bc, {t, tMin, 1}, Method -> {"EquationSimplification" -> "Residual"}]; // AbsoluteTiming 

It takes only {1148.03, Null}. Thus, we saved 9 minutes with a simple replacement t->tmax t. However, if we can write the equation in matrix form, then we can use exponential Rosenbrock-type integrators. Then we can reduce the time to 1 minute.