(*Initial Conditions*)
m = 0.050;
M = 2;
r = 0.25;
g = 9.8;
\[Epsilon] = 0.1;
\[Theta]0 = \[Pi]/3;
tf = 100;
IR = m r^2;
(*Code*)
nPend = 3;
p[i_][t_] := \[Theta][i]''[
t] + ((m r g)/
IR) Sin[\[Theta][i][t]] + \[Epsilon] ((\[Theta][i][t]/\[Theta]0)^2 -
1) \[Theta][i]'[t] + ((r m Cos[\[Theta][i][t]])/IR) y''[t]
y[t_] := -(m/(M + nPend m)) r (
Total[ Table[Sin[\[Theta][j][t]], {j, 1, nPend}]])
eqns = Thread[Table[p[i][t], {i, 1, nPend}] == 0]; (*Equations of motion*)
initpos = Thread[Table[\[Theta][i][t], {i, 1, nPend}] == 1]; (*Initial position*)
initvel = Thread[Table[\[Theta][i]'[t], {i, 1, nPend}] == 0]; (*Initial velocity*)
alleqns = Join[eqns, initpos, initvel];
slist = Table[\[Theta][i], {i, 1, nPend}]; (*What NDSolve solves for*)
NDSolve[alleqns, slist, {t, 1, tf}]