# How many ODEs can NDSolve handle?

I am trying to simulate a coupled system of N blocks and springs. I am using NDSolve. For N = 25 everything goes smoothly, however for N > 25 it doesn't work. No error message, no warning, just the usual beep.

Is there any limitation on the number of ODEs NDSolve can handle? (I'm using V9.0.1.0, Home Edition.)

UPDATE: The problem is not related with some kind of restriction on NDSolve, I just checked it with a more simpler system and n=500; NDSolve worked fine. What follows is the code; as I said before for n=25 it works well, in less than a second, however for n>26 It doesn't. Why the beep? says : The kernel Local has quit (exited) during the course of an evaluation

n = 26;
U0 = 20;
Nt = 3;
b = 0.5;
umax = 2;
tr = 40.0;

S4 = Join[{u''[t] ==
1/Nt (U0*(1 - Exp[-(t/tr)]) - u[t]) -
2*u[t]*
HeavisideTheta[2*umax - Abs[2*u[t]]] + (u[n][t] -
u[t]) + (u[t] - u[t]) - b*u'[t],
u''[t] ==
1/Nt (U0*(1 - Exp[-(t/tr)]) - u[t]) -
2*u[t]*
HeavisideTheta[2*umax - Abs[2*u[t]]] + (u[2 - 1][t] -
u[t]) + (u[2 + 1][t] - u[t]) - b*u'[t],
u''[t] ==
1/Nt (U0*(1 - Exp[-(t/tr)]) - u[t]) + (u[3 - 1][t] -
u[t]) + (u[3 + 1][t] - u[t]) - b*u'[t]
}, Table[
u[i]''[t] ==
1/Nt (U0*(1 - Exp[-(t/tr)]) - u[i][t]) -
2*u[i][t]*
HeavisideTheta[2*umax - Abs[2*u[i][t]]] + (u[i - 1][t] -
u[i][t]) + (u[i + 1][t] - u[i][t]) - b*u[i]'[t], {i, 4,
n - 1}],
{u[n]''[t] ==
1/Nt (U0*(1 - Exp[-(t/tr)]) - u[n][t]) -
2*u[n][t]*
HeavisideTheta[2*umax - Abs[2*u[n][t]]] + (u[t] -
u[n][t]) + (u[n - 1][t] - u[n][t]) - b*u[n]'[t]},
Table[u[i] == 0, {i, 1, n}],
Table[u[i]' == 0, {i, 1, n}]];
Sol = NDSolve[S4, Table[u[i], {i, 1, n}], {t, 0, 70},
Method -> "ExplicitRungeKutta"];
Plot[Evaluate[Table[u[i][t], {i, 1, n}] /. Sol], {t, 0, 70}]

• You may be running out of memory. – bbgodfrey Dec 31 '15 at 10:03
• Please provide a minimal example that can reproduce the issue. – xzczd Dec 31 '15 at 10:14
• For example: n=10; func = ToExpression@Table["f" <> ToString[i] <> "[x]", {i, 0, n}];  NDSolve[Flatten@{D[#[], x] == #[] & /@ Transpose@{func, RotateLeft@func}, (# /. x -> 0) == 1 & /@ func}, func, {x, 0, 1}] // Short  works fine for any number of ODE's (I checked up to $n=1000$). This system of ODE's is $f'_{i-1}(x)=f_i(x)$, with $f_i(0)=1$. – AccidentalFourierTransform Dec 31 '15 at 16:57
• Have you tried the menu command, Help -> Why the Beep?...? – Michael E2 Jan 1 '16 at 1:36
• @AccidentalFourierTransform You're right, I just checked it with a more simpler coupled system of ODEs and NDSolve worked fine for N=500. So the problem is not the number of equations. – Jorge Jan 1 '16 at 7:49

(1 + Tanh[p (1 - (x/xmax))])/(1 + Tanh[p])