I am trying to program a $n$-body problem simulation. To calculate the position and velocity after a time-step I want to split it in multiple 2 body problems. Now I am stuck, trying to find the velocity and position of two bodies (ignoring all the others) after one time-step with given mass, start position and velocity.

@Leon why $M$ is equal to the sum of the two masses? The law of gravitational attraction is $ F_{12}(t) = G M_1 M_2 \frac {\vec {r_{12}(t)}}{|\vec {r_{12}(t)}|^3} $
–
SofiaFeb 12 at 13:48

If your time step is small enough, you can loop over all bodies and calculate the total acceleration on each of them (component-wise). Multiply with the time step to get the change of velocity and from that you get the change of location. That way, you implicitly have covered all 2-body problems.
–
AziraphaleFeb 12 at 13:48

@Sofia This is what I got after transforming $\ddot{\vec{r_1}} - \ddot{\vec{r_2}} = \frac{Gm_2\vec{r_{12}}}{|\vec{r_{12}}|} - \frac{Gm_1\vec{r_{21}}}{|\vec{r_{12}}|}$.
–
LeonFeb 12 at 14:08

@Sofia both $\vec{r_1}$ and $\vec{r_2}$ are the location vectors of the masses. Thank you for pointing out, that this was not clear.
–
LeonFeb 12 at 15:25

Where $r_2$ is the initial distance and $r_1$ the final (so if you yould have the moon falling on the earth $r_2$ was the distance from center to center and $r_1$ the radius of the earth plus the radius of the moon.

If you want the velocities at a specific distance in free fall just replace your $r_1$ with the distance you need.

Those are the analytical solutions. If you want to solve it differentially you need to go numeric. In the following example you have 2 point masses separated by a distance of 6371 km, one the mass of the earth and the other the mass of the moon. They will meet in 889 sec and the heavy mass will move x = 77 km in the right, and the light one R-x in the left direction. If you want to solve for velocity or acceleration at time t just replace x1[t] by x2'[t] or x1''[t]: