Hello community,

I am a new user of Mathematica. I am trying to NDsolveĀ a complex heterogeneous ODE.

I want to solve for all functions (for all the n and k {n,1,nmax}, {k,1,kmax}) S$nk in terms of time, as listed at "vars" below

the hardcore commands are:

eqns = Table[{D[S$nk[n][k][t], t] == (1 - S$nk[n][k][t])*SUS$nk[n][k][t], S$nk[n][k][0] == Initiators}, {n, 1, nmax}, {k,1, kmax}];

vars = Table[S$nk[n][k][t], {n, 1, nmax}, {k, 1, kmax}];

sol = NDSolve[eqns, vars, {t, tmin, tmax}];

where SUS$nk is a function of all S$nk :

SUS$nk[n_][k_][t_] := Sum[Binomial[k, m]*(V[t])^m*(1 - V[t])^(k - m), {m, n, k}];

V[t_] := Sum[Sum[kk/z p$nnkk*S$nk[nn][kk][t], {nn, 1, nmax}], {kk, 1, kmax}];

while the Initial conditions and some extra definitions are:

z = 4;

nmax = 5;

kmax = 5;

Initiators = 0.01;

tmin = 0;

tmax = 40;

P$k = (z^k Exp[-x])/k!;

P$kk = (z^kk Exp[-x])/kk!;

\[CurlyPhi]0 = 0.18;

\[CurlyPhi]standard = \[CurlyPhi]0/100;

P$\[CurlyPhi] = 1/(\[CurlyPhi]standard Sqrt[2 \[Pi]])*Exp[-((\[CurlyPhi] - \[CurlyPhi]0)^2/(2*\[CurlyPhi]standard^2))];

p$nk = Integrate[P$\[CurlyPhi], {\[CurlyPhi], (n - 1)/k, n/k}];

p$nnkk = Integrate[P$\[CurlyPhi], {\[CurlyPhi], (nn - 1)/kk, nn/kk}];

However, I am always getting this message:

NDSolve::dsfun: "{S$nk[1][1][t],S$nk[1][2][t],S$nk[1][3][t],S$nk[1][4][t],S$nk[1][5][t]} cannot be used as a function"

Any suggestion why this is happening? I have spent a lot of time on this, and also asked colleagues, but still...