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