It is neat to see some of the old work done in the field, this looks like a pretty classic treatment of the topic. It looks like they were using a fourth-order Runge-Kutta integrator, which would likely limit long term integrations accuracy (though looks sufficient for their use case). Many algorithms I have seen typically use much higher order integration methods to beat down the accumulation of numerical error.

Source: Working on my PhD in orbital mechanics of asteroids/comets, here are my open source (python/rust) orbital integration tools: https://github.com/dahlend/kete