The foregoing method is not self-starting, so a separate initial routine is used to start the process. From the initial conditions , the respective Taylor series derivatives are formed by successive differentiation of equation (5.15), analytically. Having calculated the force and derivatives for each particle, the divided differences to third order and the time-steps are determined, and the main program is started with these quantities as if they had been obtained by the main program itself.

Since each particle is assigned its own time-step, the integration cycle itself begins
by determining the next particle *i* to be advanced, which is the one with the minimum
time-step. The coordinates and velocities of this particle are then determined to
order by using equations (5.16) and (5.17), whereupon
the current force may be obtained by direct summation. At this stage, the four past time
epochs are updated to be consistent with the definition that denotes the most
recent force evaluation. New differences are now formed, including . Together
with the new , these correction terms are combined to improve the
current coordinates and velocities to highest order, whereupon the primary coordinates
are initialized by .

New time-steps are assigned initially for all particles and at the end of each
integration cycle for particle *i*. The cycle is then repeated on the particle which now
has the minimum time-step.

**Figure:** Action trajectories for 8 galaxies in the LG and the LN (solid
lines), including N-body orbits (dotted lines). Cf. Figure 5.1.

In order to be able to compare the N-body orbits with the the AVP solution, the physical
coordinates of the `NBODY1` trajectories needs to be transformed into comoving
coordinates. Since each positions in `NBODY1` is given at time *t*, the expansion
parameter at that time has to be determined. This was done by solving
(3.33) for each time *t*, using the Runge-Kutta method. The code for this
subroutine is given in Appendix B.

Mon Jul 5 02:59:28 MET DST 1999