1. (*Initial Conditions*)
  2.  
  3. m = 0.050;
  4. M = 2;
  5. r = 0.25;
  6. g = 9.8;
  7. \[Epsilon] = 0.1;
  8. \[Theta]0 = \[Pi]/3;
  9. tf = 100;
  10.  
  11. IR = m r^2;
  12.  
  13. (*Code*)
  14.  
  15. nPend = 3;
  16.  
  17. p[i_][t_] := \[Theta][i]''[
  18. t] + ((m r g)/
  19. IR) Sin[\[Theta][i][t]] + \[Epsilon] ((\[Theta][i][t]/\[Theta]0)^2 -
  20. 1) \[Theta][i]'[t] + ((r m Cos[\[Theta][i][t]])/IR) y''[t]
  21. y[t_] := -(m/(M + nPend m)) r (
  22. Total[ Table[Sin[\[Theta][j][t]], {j, 1, nPend}]])
  23.  
  24. eqns = Thread[Table[p[i][t], {i, 1, nPend}] == 0]; (*Equations of motion*)
  25. initpos = Thread[Table[\[Theta][i][t], {i, 1, nPend}] == 1]; (*Initial position*)
  26. initvel = Thread[Table[\[Theta][i]'[t], {i, 1, nPend}] == 0]; (*Initial velocity*)
  27.  
  28. alleqns = Join[eqns, initpos, initvel];
  29.  
  30. slist = Table[\[Theta][i], {i, 1, nPend}]; (*What NDSolve solves for*)
  31.  
  32. NDSolve[alleqns, slist, {t, 1, tf}]