I have a bridge, which I want to discretize on $n$ stiff rods. To find a solution I wrote the following code in mathematica:

    n = 4;

variables = Table[ToExpression["φ" <> ToString[i]], {i, 1, n}];

xComponents = Table[Cos[variables[[i]][t]], {i, 1, n}];

yComponents = Table[-Sin[variables[[i]][t]], {i, 1, n}];

endsOfRods = Table[{xComponents[[i]], yComponents[[i]]}, {i, 1, n - 1}];

xTotal = Table[
   l (Sum[endsOfRods[[i, 1]], {i, 1, j - 1, 1}] + 1/2 endsOfRods[[j, 1]]), {j,
     2, n - 1, 1}];

yTotal = Table[
   l (Sum[If[i <= n/2, endsOfRods[[i, 2]], (-1)*endsOfRods[[i, 2]]], {i, 1, 
        j - 1, 1}] + 
      If[j > n/2, (-1)*1/2 endsOfRods[[j, 2]], 1/2 endsOfRods[[j, 2]]]), {j, 
    2, n - 1, 1}];

centreOfMass = Table[{xTotal[[i]], yTotal[[i]]}, {i, 1, Length[xTotal], 1}];

velocityCentreOfMass = 
  Table[D[centreOfMass[[i]], t].D[centreOfMass[[i]], t], {i, 1, 
    Length[xTotal], 1}];

kineticEnergy = 
  Sum[If[i == 1 || i == n, 1/2*(1/3 m l^2)*D[variables[[i]][t], t]^2, 
    1/2 m velocityCentreOfMass[[i - 1]] + 
     1/2*(1/12 m l^2)*D[variables[[i]][t], t]^2], {i, 1, n, 1}];

potentialEnergy = 
  Sum[If[i == 1 || i == n, -m g l 1/2 Sin[variables[[i]][t]] , 
     m g yTotal[[i - 1]]], {i, 1, n, 1}](*gravity*)+ 
   Sum[If[i < n/2, 1/2 k (variables[[i]][t] - variables[[i + 1]][t])^2, 
     If[i == n/2, 1/2 k (variables[[i]][t] + variables[[i + 1]][t])^2, 
      1/2 k (variables[[i + 1]][t] - variables[[i]][t])^2]], {i, 1, n - 1, 
     1}](*rotational springs*);

lagrange = kineticEnergy - potentialEnergy;

Nonconservative part

F = {0, d D[(-1)*(yTotal[[n/2 - 1]] - l/2 Sin[variables[[n/2]][t]]), t]};

r = {(xTotal[[n/2 - 1]] + l/2 Cos[variables[[n/2]][t]]), (yTotal[[n/2 - 1]] - 
     l/2 Sin[variables[[n/2]][t]])};

generalizedForces = Table[F.D[r, variables[[i]][t]], {i, 1, n/2, 1}];

Numerical solution

equations = 
  Table[If[i <= n/2, 
    D[D[lagrange, D[variables[[i]][t], t]], t] - 
      D[lagrange, variables[[i]][t]] - generalizedForces[[i]] == 0, 
    D[D[lagrange, D[variables[[i]][t], t]], t] - 
      D[lagrange, variables[[i]][t]] == 0], {i, 1, n - 1, 1}];

AppendTo[equations, (-1)*(yTotal[[n/2]] + l/2 Sin[variables[[n/2 + 1]][t]] + 
      l Sin[variables[[n]][t]]) == 0];

m = 250(*kg*);
l = 0.75(*m*);
d = 4000(*Ns/m*);
EE = 2.1*10^(11)(*Pa*);
II = 1.5*10^(-6)(*m^4*);
g = 9.81(*m/s^2*);
k = EE II/l;

AppendTo[equations, Table[variables[[i]][t] == 0 /. t -> 0, {i, 1, n, 1}]];
AppendTo[equations, 
  Table[D[variables[[i]][t], t] == 0 /. t -> 0, {i, 1, n, 1}]];

solution = 
  NDSolve[equations, Table[variables[[i]][t], {i, 1, n, 1}], {t, 0, 5}];

Now the problem with this code is, that it works perfectly in case of disretization to $n=4$ parts. But changing $n$ to something bigger ($n=10$ for example) ends rather quickly where NDSolve[] only returns the same equations that were given as input.

Are there any other methods or how should one solve a system like this?

PS: This is an upgrade to THIS PROBLEM.