The second (worse) version of your algorithm resembles the explicit
Euler scheme. It is consistent to first order only,
and it does not preserve invariants of motion such as the total energy.
In fact, the errors of the Euler scheme act as excitations,
so the orbits of your planets would spiral outward over time.
Runge-Kutta methods are designed to be consistent with the
system of differential equations to higher orders such as $4$.
Yet they are not designed to preserve the underlying Hamiltonian structure,
in particular they will cause the orbits of your planets to experience energy drift
and therefore slowly spiral inward or outward over time.
For simulations over longer time scales, you may want to avoid that.
In fact, your first integrator version does avoid energy drift.
Its explicit step of updating the velocity from (the field at) the previous position
is balanced by the formally implicit step of using the already-updated velocity
to update the position.
This is essentially an implementation of the leapfrog method,
apart from the fact that the velocities thus calculated are interpreted
more accurately as half-time samples, i.e. associated with
time $\frac{t_i+t_{i+1}}{2}$ rather than with time $t_i$.
The Verlet integrator mentioned in wltrup's answer is essentially
your leapfrog method with the velocities eliminated. To do that, the Verlet method needs
to keep track of the previous two positions instead of just one position.
(Verification is left as an exercise to the reader.)
The leapfrog and Verlet methods are symplectic integrators.
These are designed to adhere to the Hamiltonian principles underlying the system dynamics, at least for problems with separable Hamiltonians, i.e. those that can be written as
$$H(q,p) = T(p) + V(q)$$
where $q$ comprises all position coordinates, $p$ comprises the conjugate momenta,
$T$ represents kinetic energy and $V$ potential energy.
Your problem does feature a separable Hamiltonian,
so you may want to stick with symplectic integrators.
That choice rules out the use of ordinary Runge-Kutta methods.
You may however want to to increase the order of consistency of your symplectic integrator.
Then slight reductions of the time step would result in even larger gains in accuracy;
or you could afford larger time steps for the same accuracy.
The Wikipedia entry on symplectic integrators lists higher-order symplectic integration methods found by Ruth and others (1,
2, 3).
Those come with special (given, constant) coefficients
which I will label $a_1,\ldots,a_n$ and $b_1,\ldots,b_n$.
With these, one step for evolving the system from time $t$ to time $t+\delta t$
takes the form
for $i=1,\ldots,n$:
- Compute vector $F$ of all forces, given position vector $q$
- Change vector $p$ of all momenta by $b_i\,F\,\delta t$
- Compute vector $v$ of all velocities, given momentum vector $p$
- Change vector $q$ of all positions by $a_i\,v\,\delta t$
Your problem has point masses, so we can reasonably normalize each force and momentum
by the mass of the point it is applied to. Thus forces turn into accelerations,
momenta turn into velocities, and step 3 becomes a no-op.
An integration step now looks like the following:
for $i=1,\ldots,n$:
- Compute vector $g$ of all accelerations, given position vector $q$
- Change vector $v$ of all velocities by $b_i\,g\,\delta t$
- Change vector $q$ of all positions by $a_i\,v\,\delta t$
The wisdom now lies in the coefficients.
Your first-mentioned integrator has $n=1, a_1=1, b_1=1$.
For an order-$4$ symplectic scheme, Ruth has found $n=4$ and
$$\begin{bmatrix}
a_1 & a_2 & a_3 & a_4 \\
b_1 & b_2 & b_3 & b_4
\end{bmatrix}
= \frac{1}{2-\sqrt[3]{2}}
\begin{bmatrix}
\frac{1}{2} & \frac{1-\sqrt[3]{2}}{2} & \frac{1-\sqrt[3]{2}}{2} & \frac{1}{2} \\
0 & 1 & -\sqrt[3]{2} & 1
\end{bmatrix}$$
If you compare that with the referenced Wikipedia article and wonder about the reversed numbering, remember that sequencing of operator application works from right to left.
My numbering agrees with reference (3).
Note that $b_1=0$ which has the consequence that the first update of $p$ resp. $v$ is a no-op.
Therefore, each full integration step effectively begins and ends with an update to $q$,
and the values of $q$ outside the loop are not used for updating $p$ resp. $v$.
If you are not interested in those particular values of $q$, you might as well
lump $a_1$ and $a_4$ together and condense the full integration step to only three partial steps. The sequence of intermediate points would be the same, but the $q$ values outside the loop would be slightly time-shifted and give a less accurate approximation to the true position at time $t$, just as the leapfrog method has its velocities slightly time-shifted. With very coarse time steps, the noticeable effect would be that orbits would become slanted in phase space.
Update: I like to identify the leapfrog method with $(a_1, b_1) = (1, 1)$. As pointed out,
that scheme produces interleaved time samples of positions $q$ and velocities $v$. To produce time-synchronized $(q, v)$ samples, one of the intermediate steps should be split in two half-steps, and the loop body rotated so that one of the two half steps occurs at the beginning and the other at the end. This is equivalent to the scheme $(a_1, a_2; b_1, b_2) = (\frac{1}{2},\frac{1}{2};0,1)$ or $(a_1, a_2; b_1, b_2) = (1,0;\frac{1}{2},\frac{1}{2})$, which Candy/Rozmus call leapfrog and pseudo-leapfrog, respectively.
Though I'd rather emphasize the underlying equivalences, it is probably in order to follow LutzL's suggestion and name the $(a_1, b_1) = (1, 1)$ method symplectic Euler instead.
Best Answer
For each of the $N$ bodies, define its position as $x_i$ and its velocity as $v_i$, then write the differential equation for $x_i$ and $v_i$ as follows
$ \dfrac{d}{dt} x_i = v_i $
$ \dfrac{d}{dt} v_i = a_i $
Now $a_i = \dfrac{1}{m_i} F_i $
where $F_i = \displaystyle \sum_{ j , j \ne i} \dfrac{ G m_i m_j (x_j - x_i) }{\| x_i - x_j \|^3} $
The state vector of the whole system is defined as
$ x = [x_1, x_2, x_3, ..., x_N , v_1, v_2, v_3, ..., v_N ]^T$
Then you have
$ \dot{x} = f(x) $
And you just have to integrate this vector differential equation numerically, using an accurate method like the Runge-Kutta $4$-th order method. This will give you the position of all the $N$ bodies as a function of time.