N-Body problem (Taylor series)

[Comp.parallel moderator:
This is a bit of a non-standard cross post, but it falls within the
guidelines of c.p. I am going to set up the Followup line to other
unmoderated groups. People willing to post a followup should consider
further editing the header. If someone really wants c.p. in the answer,
I will consider keeping the full original Newsgroup line.
It just has to be c.p. relevant.
--the c.p. moderator ]


We're writing project solving N-Body problem. The method we're using
for aproximation of motion equations is based on Taylor series.
* To count derivates (like accelaration and so on) we use scheme of
Dybczynski (Acta Astronomica, 1979); by the way, have you ever heard
about this man?
* To set length of integration step (time interval) we've implemented
adaptive Sitarski's method, which is based on comparing two last
derivates of ONE body. How to choose THAT one body? Should I choose
maximum force from Keppler's equotation?
But... we meet some problems, which were not supposed to appear in
this method:
The coordinates (x,y,z) and velocities are computed independly and
meet something like slingshot effect but it should not appear in this
method (like author of the method was kind to said us). The Mars comes
near by Sun and doesn't manage to slow down.
Afterwards I've tried to turn off the gravity influence for all bodies
if the step length is less than 0.0001 AU. When step length is
supposed to be bigger I turn it on again.
We've check all input data. Settings we tested includes: accuracy-
epsilon = 10^-14, order of Taylor series = 20.
Could any one help us? Any ideas what's going on wrong?

Maciej.
(dzieciou)