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