All multibody orbit simulations will become inaccurate no matter what the method of calculation if you integrate it far enough into the future. This is due to the accumulation of round-off error. People have known about this for over a centry. The best numerical method for integrating the equations of motion directly (Cowell's method) is the Gauss-Jackson predictor-correction. However even with the Gauss-Jackson predictor-correction, the round off error accumulates fairly rapidly for Cowell's method of orbit integration. If the forces perturbing the trajectory are small enough then the best approach is to use a "perturbation method". A fourth/fifth order Fehlberg Runge-Kutta integrator is arguably the best approach for integrating a perturbation method (that's what I use). The Fehlberg Runge-Kutta is a hybrid method that uses two different Runge-Kutta schemes to calculate truncation error and thus determine the optimal step size. Many books have been written about perturbation methods for orbital mechanics.

Gary