First, good article, I was wondering many of this myself. I had a couple thoughts on reading it, from my experience in doing numerical integration. 1. With regards to integration schemes, given how smooth planetary orbits are, you could probably use a much higher order integration scheme. The RK45 that you used is good, I remember using a 7th/8th order Runga-Kutta method that I liked too. Or if you wanted to be super serious, you could use ODEPACK, a standard numerical library for solving ODEs. I remember being impressed on the speed gain I would get using ODEPACK compared with a simpler Runga-Kutta. Whether it is worth the programing time is another question. 2. As for not being able to find stable Lagrangian Points, I would want to see more simulation results, but I would guess that the numerical errors from your integration scheme are causing them to loose stability. Runga-Kutta methods do not conserve energy, so for physical systems where that is important (like planetary dynamics), they will introduce errors in the long term. These effects are often worst around stable points. There exist integration schemes that do preserve energy, but these are generally much more complex. Other possibilities are that perturbations from Minimus are causing you problems, or possibly that you have a bug somewhere (which I'm sure isn't the case, after all, my programs _never_ have bugs, why would anybody else's? :-)) p.s. This is the first time I've posted here!