Use Mathematica to investigate the limit cycles, if they exist, of the system

˙ x = y − ǫ(a1x + a2x2 +· · ·+aMxM), ˙ y = −x, as the parameter ǫ varies from zero to infinity. When ǫ is large the step size has to be small; in fact, it would be better to use a variable step length in the numerical solver.