The "leap frog" method is very efficient and stable, and works well for any dynamic particle/field system, including plasmas. For gravity, it's simple. Here is all you do for a single iteration on a single planet (1-body problem, single planet around stationary sun):
public void Push()
{
Position += Gravity.dT * Velocity;
Velocity += Gravity.dT * GravityAccelerationVector(Position);
}
Where "Gravity.dT" is a uniform time step, in an arbitrary measure of time. I'm using System.Windows.Vector, but any custom Vector class will do, as long as it supports basic multiplication and addition. The trick is that the Position and Velocity aren't at the "same time", which is very common for most integration methods. Rather, they're staggered. The Position on iteration N is updated based on the Velocity from iteration N - 1/2, but then the Velocity at iteration N+1/2 is updated based on the Position at iteration N.
The N-body version looks like this:
public static void PushPlanets(Planet[] planets)
{
// Position Push at iteration N + 0:
foreach(var p in planets)
p.Position += Gravity.dT * p.Velocity; // Velocity from N - 1/2
// Velocity Push at iteration N + 1/2:
foreach (var p in planets)
{
Vector TotalGravity = new Vector(0,0);
foreach (var pN in planets)
{
if (pN == p) continue;
TotalGravity += pN.Mass * p.Mass * GravityAccelerationVector(p.Position - pN.Position);
}
TotalGravity += Sun.Mass * p.Mass * GravityAccelerationVector(p.Position); // Solar acceleration
p.Velocity += Gravity.dT * TotalGravity;
}
Where
public static Vector GravityAccelerationVector(Vector position)
{
return Vector.Multiply(-G / position.LengthSquared / position.Length, position);
}
The N-body is only more complicated because instead of a single gravitational source, there are several. The code format is the same, though: each planet's position gets pushed by the N-1/2 Velocity, then we figure out the total gravity acceleration on each planet based on the new positions, then we push each planet's Velocity by that total acceleration.
Other methods, even high order methods, are often unstable because they linearly project the next step based on position and velocity at the same time. This will always err on the side of adding energy to the system, and the orbits will gradually move further and further outward. Other methods can overcompensate for this natural error and remove energy from the system. In general, one wants an energy neutral solution. The leap frog method will gradually err in terms of phase of the orbits, but not in overall energy.