Euler integration is the simplest possible integration of motion equations and after implementing it you will find that small values are allowed. Problems with stiffness occurs too and it is well known disadvantage of low order schemes of ODE integration. To make solution more accurate and stable you should consider using more complex integrator. It may be an explicit scheme from Runge-Kutta family (second order Mid-Point method could be good choice) or one of unconditionaly stable implicit schemes (i.e. Backward Euler). Besides I am not a fan of implicit schemes (complex solution + huge computational cost) I propose to use Heun Predictor-Corrector scheme. It is rather easy to implement (if you understand Euler integration) and easy to understand (it is also important to know how your integrator works, not only how your integrator is implemented).

We consider problem of ODE:

(12) |

Detailed description and discussion about semi-implicit schemes you can find in great book by M. G. Ancona [2]. Heun's predictor/corrector scheme consist of two main steps, a predictor:

(13) |

and corrector step:

Whole procedure is quite easy - first we compute forces, then a simply Euler step is done (prediction). Result of that is substituted into corrector as an argument of force computation procedure. Then we use correction central scheme to update and finalize integration. It is easy, cost us only one additional step of force computation and gives semi implicitness which result in second order accuracy of the solution.