Numerics for the spin orbit equation of Makarov with constant   eccentricity
  We present an algorithm for the rapid numerical integration of a time-periodic ODE with a small dissipation term that is $C^1$ in the velocity. Such an ODE arises as a model of spin-orbit coupling in a star/planet system, and the motivation for devising a fast algorithm for its solution comes from the desire to estimate probability of capture in various solutions, via Monte Carlo simulation: the integration times are very long, since we are interested in phenomena occurring on times similar to the formation time of the planets. The proposed algorithm is based on the High-order Euler Method (HEM) which was described in~\cite{hem}, and it requires computer algebra to set up the code for its implementation. The pay-off is an overall increase in speed by a factor of about $7.5$ compared to standard numerical methods. Means for accelerating the purely numerical computation are also discussed 
Michele V. Bartuccelli, Jonathan H.B. Deane and Guido Gentile
math.CA