Euler-Richardson or midpoint method is good for constant acceleration.

Normally, I use Runge-Kutta 4th order method, it is rough 100 times better resolution than midpoint error. This is good for most cases.

Sometime, I used Runge-Kutta Fehlberg adjustable time step method for interaction for like 1/r[sup]2[/sup] (when r is small , the error become larger).

However, you need to asjust time step carefully.

Smaller time will give you smaller run time error until run-offer error become major problem (10[sup]-6[/sup] for float, 10[sup]-13[/sup] for double).

If the particle size is very small , or there are many particels. the error will accumuate.

You need to study numerical method to understand what is the best way to used for different cases.

I would suggest you use Runge-Kutta 4th order method for most of the cases.

The error could be due to numerical method. However, it could be due to the coding.

You can calculate error if you know the theoretical value in advance.

Otherwise, it will depend on how much you know about the system to have a better guess.