[1] http://www.damtp.cam.ac.uk/user/tong/relativity/dynrel.pdf
That's an apples-to-oranges comparison, though. The 10^-70 number is a relative error while the planck length is just a length. The numerical computations in the paper were likely done in a system of "simulation units" where the actual lengths, velocities, forces, etc. are normalized to values of order unity. This has advantages for the numerical aspect, in terms of preserving floating point accuracy. But it also means that to translate it to a physical system (e.g., a triple system of stars), the simulation needs to be scaled to physical units. The 10^-70 that's quoted just means that the values should be numerically accurate to 1 part in 10^70.
If you wanted to compare the (floating point) accuracy of the simulation with the Planck length, you would need to scale the simulation to a physical size and see what the 10^-70 fractional error would translate to in physical units.
> At first, using the obtained initial conditions, we checked the 137 periodic orbits by means of the high-order Taylor series method in the 100-digit precision with truncation errors less than 10^−70 , and guaranteed that they are indeed periodic orbits.
> Besides, we use the CNS with even smaller round-off error (in 120-digit precision) and truncation error (less than 10^−90 ) to guarantee the reliability of these 27 families.
They do reference some software (e.g., "dop853"), but I'm not familiar with the details of those ODE solvers.