Efficiently tabulating system of coupled ODEs for NDSolve

I am solving a (large) system of second order coupled nonlinear ODEs. I have been increasing the number of dependent variables present in the system slowly, to converge on more physically realistic results, but have experienced some pretty severe bottlenecks in the code.

I will present some theory and my code below, but to be clear, my question is this: What is the best (i.e. most efficient) way to tabulate a large system of indexed functions, to be passed into NDSolve?

I've tried a few things to improve efficiency which I'll list after I present the code, but first some details.

The system under consideration here is a generalization of the one I discussed here and in particular here. Qualitatively, I am trying to solve a system of coupled nonlinear (up to order 4 in the dependent variables) ODEs which can take on complex values.

The governing equation for the $n-th$ ODE is of the form:

$\sum_{k=-M}^M (Q_{nk}+Q_{kn}) \ddot{a}_k -\sum_{\ell=-M}^M\sum_{k=-M}^M S_n(k,\ell) \dot{a}_{\ell} \dot{a}_k + U_n = 0$ where

$Q_{nk}=\frac{1}{4}\sum_{j=1}^MP_{jn}P_{(-j)k}$,

with

$P_{mn}=\frac{1}{|n|\sqrt{|m|}}(a_{m-n}-\frac{a_ma_{-n}}{2})$,

(and the contsraints that $a_{m-n}=0$ for $|m-n|>M$ and $a_{m-n}=0$ when $Sign(m-n)=Sign(n)$),

$S_n(k,\ell)=\left(\frac{\partial Q_{k\ell}}{\partial a_n}-\frac{\partial Q_{n\ell}}{\partial a_k}-\frac{\partial Q_{\ell n}}{\partial a_{k}}\right)$ and

$U_n = \frac{\partial V}{\partial a_n}$ where

$V = \frac{1}{8}\left(\sum_{n=1}^M\frac{a_na_{-n}}{n}\right)^2 +\frac{1}{2}\sum_{n=1}^M\frac{a_na_{-n}}{n^2} -\frac{1}{2} \left(\sum_{n=1}^M\frac{a_na_{-n}}{n}\right)\left(\sum_{n=1}^M\frac{a_na_{-n}}{n}\right) + \sum_{i+j+k=0,(i,j,k)\neq 0}\frac{a_ia_ja_k}{4ij}$

and $M$ is the highest mode under consideration.

Furthermore, I was experiencing problems using NDSolve when the functions take on complex values (as discussed here) so that I have decomposed each dependent variable as $a_n =AR_n+i*AI_n$, where $AR_n,AI_n$ are real functions. For M small I am finding coherent results with my code, however, several bottlenecks are impeding my progress to higher M. At this point it makes sense to list my code:

M = 10; (*Number of modes*)
cond1[r_, s_] :=
cond1[r, s] =
If[Sign[s] == Sign[r - s], 0, 1];
cond2[r_, s_] :=
cond2[r, s] =
If[Abs[r - s] > M, 0, 1]; F[n_] :=
F[n] = If[n == 0, 1, Abs[n]];
a[0][t_] := a[0][t] = 1;
P[m_, k_, t_] :=
P[m, k, t] =
1/(Sqrt[F[m]]*F[k])*(cond1[m, k]*cond2[m, k]*a[m - k][t] -
1./2.*a[m][t]*a[-k][t]) ;
Q[k_, l_, t_] :=
Q[k, l, t] = 1./4.*Sum[P[m, k, t]*P[-m, l, t], {m, 1, M}];
q[m_, k_, n_, t_] :=
q[m, k, n, t] =
1./(Sqrt[F[m]]*F[k])*(cond1[m, k]*cond2[m, k]* Boole[m - k == n] -
a[m][t]/2.*Boole[-k == n] - a[-k][t]/2.*Boole[m == n]);
SQ[k_, l_, n_, t_] :=
SQ[k, l, n, t] =
1./4.*Sum[
P[-m, l, t]*q[m, k, n, t] + P[m, k, t]*q[-m, l, n, t], {m, 1, M}];
SS[k_, l_, n_, t_] :=
SS[k, l, n,
t] = (SQ[k, l, n, t] - SQ[n, l, k, t] - SQ[l, n, k, t]);
Fu[a_, b_, c_] := Fu[a, b, c] = If[a + b + c == 0, 1, 0];
Fu2[a_] := Fu2[a] = If[a == 0, 0, 1];
V[t_] := V[
t] = (Sum[-a[n][t]*a[-n][t]/(2*n), {n, 1, M}]^2 +
2.*Sum[a[n][t] a[-n][t]/(4*n^2), {n, 1, M}] +
Sum[-a[n][t]*a[-n][t]/(2 n), {n, 1, M}]*
Sum[a[n][t] a[-n][t]/n, {n, 1, M}] +
Sum[Fu[n, m, o]*Fu2[n]*Fu2[m]*a[n][t]*a[m][t]*a[o][t]*
Abs[o]/(8*Abs[F[n]*F[m]*F[o]]), {n, -M, M}, {m, -M, M}, {o, -M,
M}]);
Join[Table[With[{i = i}, U2[i] = D[V[t], a[i][t]]], {i, 1, M}],
Table[With[{i = i}, U2[-i] = D[V[t], a[-i][t]]], {i, 1, M}]];
Gov[n_, t_] :=
Sum[a[l]''[t]*(Q[n, l, t] + Q[l, n, t]), {l, -M, M}] -
Sum[SS[k, l, n, t] a[k]'[t] a[l]'[t], {k, -M, M}, {l, -M, M}] +
1/2*U2[n];
Table[With[{i = i}, a[i][t_] := ar[i][t] + I*ai[i][t]], {i, 1, M}];
Table[With[{i = i}, a[-i][t_] := ar[i][t] - I*ai[i][t]], {i, 1, M}];
Table[ar[i][t] \[Element] Reals, {i, 1, M}];
Table[ai[i][t] \[Element] Reals, {i, 1, M}];
Table[
With[{i = i}, Test[i] = Evaluate[Gov[i, t]]], {i, 1, M}];
Table[With[{i = i}, ReGov[i, t_] := ComplexExpand[Re[Test[i]]]], {i,
1, M}];
Table[With[{i = i}, ImGov[i, t_] := ComplexExpand[Im[Test[i]]]], {i,
1, M}];
ao = 0.1;
eqns = Join[Parallelize[Table[ReGov[i, t] == 0, {i, 1, M}]],
Parallelize[Table[ImGov[i, t] == 0, {i, 1, M}]], Table[ar[i][0] == ao, {i, 1}],
Table[ar[i][0] == 0, {i, 2, M}], Table[ai[i][0] == 0, {i, 1, M}],
Table[ar[i]'[0] == 0, {i, 1, M}], Table[ai[i]'[0] == ao, {i, 1}],
Table[ai[i]'[0] == 0, {i, 2, M}]];
MMode = NDSolve[eqns,
Join[Table[ar[i][t], {i, 1, M}], Table[ai[i][t], {i, 1, M}]], {t,
0, 10}];


The first bottleneck occurs when defining a variable that evaluates the governing equations:

Table[
With[{i = i}, Test[i] = Evaluate[Gov[i, t]]], {i, 1, M}];


Note, for the M=10 case the Timing on my machine is about 1.55, while for the M=20 case, the Timing is about 28.555. In particular, the bottleneck seems to come from evaluating the function SS[n,k,l,t] in the definition of Gov[i,t]. There are quite a few terms in this function, but they all involve low order polynomials of maximum degree 4 in the dependent variables. I'm not sure if there is a way to speed these calculations up, or if this is a necessary area of intense computation.

As mentioned above, I decompose my governing equations into real and imaginary parts, to speed up the NDSolve portion of the computation. This leads to the next, more severe, bottle neck, which comes from tabulating the real and imaginary parts of the governing equations, i.e. defining the variable "eqns". I have used ComplexExpand[Re(Im)[*]] to get the real and imaginary parts of the governing equations, respectively, but this seems to be computationally expensive and accounts for the gross majority of the computation time. The reason for this post is to find out if there is a more efficient way to do this. Note, for the M=5 case the Timing of the eqns calculation is 0.0173 while for the M=10 case it is 12.25, i.e. 4 orders of magnitude larger!

Qualitatively, I wonder if I have not effectively constrained some of the definitions in my code, and if this generality is leading to significant slow down. Also, I feel as if I need to gain some understanding of the compile function, which seems to be very useful when defining functions, but the subtleties of using it with indexed functions to be passed into NSDSolve is beyond my very naive understanding and I cannot find any relevant examples to study. Finally, an alternative approach would be using matrices instead of term by term operations, but I'm not convinced this will lead to a significant speed up in Mathematica, but would welcome any advice on this.