# 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[1]''[t] ==
1/Nt (U0*(1 - Exp[-(t/tr)]) - u[1][t]) -
2*u[1][t]*
HeavisideTheta[2*umax - Abs[2*u[1][t]]] + (u[n][t] -
u[1][t]) + (u[2][t] - u[1][t]) - b*u[1]'[t],
u[2]''[t] ==
1/Nt (U0*(1 - Exp[-(t/tr)]) - u[2][t]) -
2*u[2][t]*
HeavisideTheta[2*umax - Abs[2*u[2][t]]] + (u[2 - 1][t] -
u[2][t]) + (u[2 + 1][t] - u[2][t]) - b*u[2]'[t],
u[3]''[t] ==
1/Nt (U0*(1 - Exp[-(t/tr)]) - u[3][t]) + (u[3 - 1][t] -
u[3][t]) + (u[3 + 1][t] - u[3][t]) - b*u[3]'[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[1][t] -
u[n][t]) + (u[n - 1][t] - u[n][t]) - b*u[n]'[t]},
Table[u[i][0] == 0, {i, 1, n}],
Table[u[i]'[0] == 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. Dec 31, 2015 at 10:03
• Please provide a minimal example that can reproduce the issue. Dec 31, 2015 at 10:14
• For example: n=10; func = ToExpression@Table["f" <> ToString[i] <> "[x]", {i, 0, n}];  NDSolve[Flatten@{D[#[[1]], x] == #[[2]] & /@ 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$. Dec 31, 2015 at 16:57
• Have you tried the menu command, Help -> Why the Beep?...? Jan 1, 2016 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. Jan 1, 2016 at 7:49

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