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

Here is a way to type the initial data:

initialData =
 {\[Theta][1, 0] == 1, Derivative[0, 1][\[Theta]][1, 0] == 2,
  \[Theta][2, 0] == -1, Derivative[0, 1][\[Theta]][2, 0] == 4,
  \[Theta][3, 0] == 3, Derivative[0, 1][\[Theta]][3, 0] == 0}

This is a difficult task for a beginner. Here is an attempt:

r[n_, t_] = \[Zeta][t] {0, 0, 1} + 
   Sum[{1, 0, 0} l[i] Sin[\[Theta][i, t]] -
     {0, 0, 1} 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]];
eqs = Block[{nn = 2, m = Function[1],
    l = Function[1], cd = 1, \[Zeta] = Function[1], g = 1},
   Table[eqLag[n], {n, nn}]] // Simplify
sols = NDSolveValue[{eqs,
   \[Theta][1, 0] == 1, \[Theta][2, 0] == 1,
   Derivative[0, 1][\[Theta]][1, 0] == 1,
   Derivative[0, 1][\[Theta]][2, 0] == 1},
  {\[Theta][1, t], \[Theta][2, t]}, {t, 0, 4}]
ParametricPlot[sols, {t, 0, 4}]