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}]