To a first approximation, this behaves just like any other mass-under-a-force problem: since the gravitational force is
F=G*m1*m2/r^2, the acceleration on body 1 is
a1=G*m2/r^2 and that on body 2 is
a2=G*m1/r^2 — note that these accellerations are in the direction of the other body.
To keep this pseudo-code relatively minimal, I'm going to presume that there's already a well-defined
Vector class with all of the relevant operations (addition, subtraction, scalar multiplication and normalization) available, and all of my vector variables will be prefaced with a 'v'; e.g.,
Given this setup, then a first approximation is pretty straightforward — but there's one huge caveat that I'll cover at the end.
First of all, there's a straightforward initialization step:
vPos1 = vInitPos1; vVel1 = vInitVel1;
vPos2 = vInitPos2; vVel2 = vInitVel2;
Then on every update, we first find the vector from the first body to the second — of course, the vector from the second body to the first is just the negative of this. While we're at it, we pre-scale it by the gravitational constant and the inverse squared distance between the two vectors; these values will be the same for both accelleration computations, so now's a good time to get them out of the way. Note that since we'll normalize vDeltaPos by dividing it by its length, and then we divide it again by its squared length, we could actually just divide by the length cubed - but vector normalization routines often have special code for numerical stability, so it's generally best to do this 'the hard way'.
vDeltaPos = vPos2-vPos1;
radSq = vDeltaPos.LengthSquared();
vScaledDeltaPos = vDeltaPos*(GRAVITATIONAL_CONSTANT/radSq);
Now with this in hand, the acceleration for the two bodies is easy; body 1's accel just points towards body two, scaled by body 2's mass, while body 2's accel points towards body 1, scaled by that body's mass:
vAccel1 = mass2*vScaledDeltaPos;
vAccel2 = (-mass1)*vScaledDeltaPos;
With the accelerations in hand, then the update step is just two iterations of the Euler method, updating velocity from acceleration and updating position from velocity:
vVel1 += deltaTime*vAccel1;
vVel2 += deltaTime*vAccel2;
vPos1 += deltaTime*vVel1;
vPos2 += deltaTime*vVel2;
And that's all there is to it! But... there's a huge caveat here. Actually, two huge caveats. The first is that Euler's method is as inaccurate as a differential equation solver can get while still being useful; there are methods that give you much better accuracy for a given delta-Time, but they're a lot more complicated. (If you're curious about these, the magic names to start with are Runge-Kutta).
The second, though, is a little more subtle, and the code above luckily accounts for it. You may notice that I update the velocity before I update the position; this means that the position update uses the updated velocity rather than the 'old' velocity. This minor tweak helps to correct for a problem with the naive Euler methods: they don't necessarily conserve the energy of the system, so that even a configuration that looks like it should be a stable mutual orbit can wind up with the two bodies flying off towards infinity. Using the updated velocity is a very simple version of a so-called semi-implicit Euler method, which helps with energy conservation (at least to the approximation provided by the timestep), but if you want guaranteed long-term conservation - i.e., stable orbits over long timescales - then you'll want to look even deeper down the rabbit hole, into methods that are guaranteed to be energy-preserving.