'Repa' performance for planetary simulation

Most Euler numeric integration methods suffer from cumulative round-off error that will eventually cause the simulation to "blow up". You may want to investigate advanced numerical integration methods, such as 4th-order Runge-Kutta or predictor-corrector.

Another place where n-body problem simulations become sticky is when two bodies get very close, such as a moon with a very eccentric orbit about its planet. If one uses fixed time increments for the simulation, the error during large changes of angular velocity can lead to division-by-zero errors or division by very small values that result in the simulation blowing up. Use of a variable delta-t that depends on the angular velocity can be beneficial.

These suggestions are based on running many such simulations as a project for an undergraduate physics course I took in 1973, while testing various numerical integration methods. Runge-Kutta and predictor corrector methods have been around since the dawn of digital computing and a number of books are available. See, e.g., Numerical Recipes: The Art of Scientific Computing by William H. Press, Brian P. Flannery, Saul A. Teukolsky and William T. Vetterling. (Cambridge University Press, 1989)