How to solve numerically this system of vector differential equations (along with its initial conditions) without using any physical trick (like the reduced mass $\mu=\frac{m_1m_2}{m_1+m_2}$):
$$G\frac{\require{cancel}\cancel{m_1}m_2}{||\vec{r}_1-\vec{r}_2||^2}\frac{\vec{r}_2-\vec{r}_1}{||\vec{r}_1-\vec{r}_2||}=\require{cancel}\cancel{m_1}\frac{d^2\vec{r}_1}{dt^2} \tag{1}$$
$$G\frac{m_1\require{cancel}\cancel{m_2}}{||\vec{r}_1-\vec{r}_2||^2}\frac{\vec{r}_1-\vec{r}_2}{||\vec{r}_1-\vec{r}_2||}=\require{cancel}\cancel{m_2}\frac{d^2\vec{r}_2}{dt^2} \tag{2}$$
where
$$\vec{r}_1=\vec{r}_1(t)=x_1(t)\hat{i}+y_1\hat{j}+z_1(t)\hat{k}$$ $$\vec{r}_2=\vec{r}_2(t)=x_2(t)\hat{i}+y_2\hat{j}+z_2(t)\hat{k}$$
This is my attempt to solve it:
- Separate into scalar equations:
$$\left\{\begin{matrix} \frac{Gm_2}{[(x_1-x_2)^2+(y_1-y_2)^2+(z_1-z_2)^2)]^{3/2}}(x_2-x_1)=\ddot{x}_1\\ \frac{Gm_2}{[(x_1-x_2)^2+(y_1-y_2)^2+(z_1-z_2)^2)]^{3/2}}(y_2-y_1)=\ddot{y}_1 \\ \frac{Gm_2}{[(x_1-x_2)^2+(y_1-y_2)^2+(z_1-z_2)^2)]^{3/2}}(z_2-z_1)=\ddot{z}_1 \end{matrix}\right.$$
$$\left\{\begin{matrix} \frac{Gm_1}{[(x_1-x_2)^2+(y_1-y_2)^2+(z_1-z_2)^2)]^{3/2}}(x_1-x_2)=\ddot{x}_2 \\ \frac{Gm_1}{[(x_1-x_2)^2+(y_1-y_2)^2+(z_1-z_2)^2)]^{3/2}}(y_1-y_2)=\ddot{y}_2 \\ \frac{Gm_1}{[(x_1-x_2)^2+(y_1-y_2)^2+(z_1-z_2)^2)]^{3/2}}(z_1-z_2)=\ddot{z}_2 \end{matrix}\right.$$
- Reduce order:
- I define the new variables $$\begin{matrix} u_1 \triangleq x_1 & u_2 \triangleq \dot{x}_{1}=\dot{u}_{1} & \dot{u}_{2}=\ddot{x}_1\\ u_3 \triangleq y_1 & u_4 \triangleq \dot{y}_{1}=\dot{u}_{3} & \dot{u}_{4}=\ddot{y}_1\\ u_5 \triangleq z_1 & u_6 \triangleq \dot{x}_{1}=\dot{u}_{5} & \dot{u}_{6}=\ddot{x}_1\\ u_7 \triangleq x_2 & u_8 \triangleq \dot{x}_{2}=\dot{u}_{7} & \dot{u}_{8}=\ddot{x}_2\\ u_9 \triangleq y_2 & u_{10} \triangleq \dot{y}_{2}=\dot{u}_{9} & \dot{u}_{10}=\ddot{y}_2\\ u_{11} \triangleq z_2 & u_{12} \triangleq \dot{z}_{2}=\dot{u}_{11} & \dot{u}_{12}=\ddot{z}_2\\ \end{matrix}$$
- I write the system of diferential equations in the form $\vec{u}'(t)=\vec{F}(\vec{u},t)$ with
$$\vec{u}'(t)=\begin{bmatrix} \dot{u}_{1}\\ \dot{u}_{2}\\ \dot{u}_{3}\\ \dot{u}_{4}\\ \dot{u}_{5}\\ \dot{u}_{6}\\ \dot{u}_{7}\\ \dot{u}_{8}\\ \dot{u}_{9}\\ \dot{u}_{10}\\ \dot{u}_{11}\\ \dot{u}_{12}\\ \end{bmatrix}$$
$$\vec{F}(\vec{u},t)=\begin{bmatrix} {u}_{2}\\ \frac{Gm_2}{[(u_1-u_7)^2+(u_3-u_9)^2+(u_5-u_{11})^2)]^{3/2}}(u_7-u_1)\\ {u}_{4}\\ \frac{Gm_2}{[(u_1-u_7)^2+(u_3-u_9)^2+(u_5-u_{11})^2)]^{3/2}}(u_9-u_3)\\ u_{6}\\ \frac{Gm_2}{[(u_1-u_7)^2+(u_3-u_9)^2+(u_5-u_{11})^2)]^{3/2}}(u_{11}-u_5)\\ u_{8}\\ \frac{Gm_1}{[(u_1-u_7)^2+(u_3-u_9)^2+(u_5-u_{11})^2)]^{3/2}}(u_1-u_7)\\ u_{10}\\ \frac{Gm_1}{[(u_1-u_7)^2+(u_3-u_9)^2+(u_5-u_{11})^2)]^{3/2}}(u_3-u_9)\\ u_{12}\\ \frac{Gm_1}{[(u_1-u_7)^2+(u_3-u_9)^2+(u_5-u_{11})^2)]^{3/2}}(u_5-u_{11})\\ \end{bmatrix}$$
- Solve the system using RK4(5), for example, in Matlab, using the function ode45, with the initial condition vector $\vec{u}_0$ and the time span $tspan$ from $t0=0$ to $tmax=31536000$ (1 year).
G=6.67408e-11; %gravitational constant m^3 kg^-1 s^-2
m1=1.98855e30; %mass sun kg
m2=5.97237e24; %mass earth kg
r0=149.6e9; %distance sun-earth m
v0=29.78e3; %earth-to-sun speedm/s
r = @(U) ((U(1)-U(7))^2+(U(3)-U(9))^2+(U(5)-U(11))^2)^(3/2);
dUUdt = @(t,U,G,m1,m2) [U(2);...
G * m2 / r(U) * (U(7)-U(1));...
U(4);...
G * m2 / r(U) * (U(9)-U(3));...
U(6);...
G * m2 / r(U) * (U(11)-U(5));...
U(8);...
G * m1 / r(U) * (U(1)-U(7));...
U(10);...
G * m1 / r(U) * (U(3)-U(9));...
U(12);...
G * m1 / r(U) * (U(5)-U(11))];
dUdt = @(t,U) dUUdt(t,U,G,m1,m2);
U0 = zeros(12,1);
U0(7) = r0; %x0
U0(10) = v0; %vy0
tspan=0:86400:31536000; % 1 year in seconds with a step of 1 day
[T,Y]=ode45(dUdt,tspan,U0);
plot(Y(:,7),Y(:,9)) % corroboration that orbit is aprox circular
(Y(:,7).^2+Y(:,9).^2).^(1/2) %corroboration that earth-sun distance is aprox constant
(Y(:,8).^2+Y(:,10).^2).^(1/2) %corroboration that earth-to-sun speed is aprox constant
It's quite simple actually. You strat by wtiting the second order dfferential equation into first order equations. Let's start by considering the x ccordinate here, equation can be re-written as
$ \vec {\ddot x} = \dfrac{1}{m}\bigg(G\dfrac{m_1m_2}{(||\vec{x_1}-\vec{x_2}||)^2}-\dfrac{\vec{x_2}-\vec{x_1}}{||\vec{x_1}-\vec{x_2}||} \bigg)$
Let $\dot x= v_x$ and $\ddot x = \dot v_x$
which will give you set of two equations written as
$\dot x= v_x$
and
$\dot v_x = \dfrac{1}{m}\bigg(G\dfrac{m_1m_2}{(||\vec{x_1}-\vec{x_2}||)^2}-\dfrac{\vec{x_2}-\vec{x_1}}{||\vec{x_1}-\vec{x_2}||} \bigg)$
and your initial conditions become
$x(0) = x_0\\ \dot v_x(0) = \dot x_0$
same thing can be done for $y$ and $z$ coordinates for each equation, and then you can use any forward integration scheme such as Euler's or Runge-Kutta to integrate numerically till desired time.
Hope this is helpful enough