English: Done in
Mathematica with the following code:
<< VariationalMethods`
tmax = 1000;
solution = Block[{m, g, \[Theta], L, x, y, v, t},
x = Accumulate[
Sin[{Subscript[\[Theta], 1], Subscript[\[Theta], 2]}]];
y = Accumulate[ -Cos[{Subscript[\[Theta], 1], Subscript[\[Theta],
2]}]];
v = D[{x, y} /.
Subscript[\[Theta], i_] :> Subscript[\[Theta], i][t], t];
L = Plus @@ (1/
2 Map[# . # &, Transpose[v]] - (y /.
Subscript[\[Theta], i_] :> Subscript[\[Theta], i][t]));
NDSolve[Join[Simplify[
Table[
EulerEquations[L, Subscript[\[Theta], i][t], t], {i, 1, 2}
]
], {Subscript[\[Theta], 1][0] == 22/7,
Subscript[\[Theta], 2][0] == 22/7,
Subscript[\[Theta], 1]'[0] == Subscript[\[Theta], 2]'[0] ==
0}], {Subscript[\[Theta], 1], Subscript[\[Theta], 2]}, {t, 0,
tmax}]
];
ParametricPlot[
Table[Subscript[\[Theta], i][t], {i, 2}] /. solution, {t, 0, tmax}]