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][1]. 


  [1]: https://mathematica.stackexchange.com/questions/110495/how-to-systematically-approach-when-solving-a-large-system-of-differential-equat?noredirect=1#comment299753_110495