Here is a way:
r[n_, t_] = \[Zeta][t] UnitVector[3, 3] +
Sum[UnitVector[3, 1] l[i] Sin[\[Theta][i, t]] -
UnitVector[3, 3] l[i] Cos[\[Theta][i, t]], {i, n}];
rdot[n_, t_] = D[r[n, t], t];
tt = 1/2 Sum[m[n] Total[rdot[n, t]^2], {n, nn}];
vv = -g*Sum[
m[n] Sqrt[Total[r[n, t]^2]] Cos[\[Theta][n, t]], {n, nn}];
ll = tt - vv;
rr = 1/2 cd*Sum[Total[rdot[n, t]^2], {n, nn}];
eqLag[k_] =
D[D[ll, D[\[Theta][k, t], t]], t] -
D[ll, \[Theta][k, t]] == -D[rr, D[\[Theta][k, t], t]];
Block[{nn = 3, m = Function[1], l = Function[1],
cd = 1, \[Zeta] = Function[1], g = 1},
eqs = Table[eqLag[n], {n, nn}] // Simplify;
initialData =
Table[D[\[Theta][k, t], {t, order}] == 1, {k, nn}, {order, 0,
1}] /. t -> 0 // Flatten;
sols = NDSolveValue[Flatten@{eqs, initialData},
Table[\[Theta][k, t], {k, nn}],
{t, 0, 4},
Method -> {"EquationSimplification" -> "Residual"}]]
