In order to have the three-body in 3D, the z-coordinate has to be implemented. Using the equations from the three-body we added the z-cooridate to all the equations. This resulted in a system of 18 ODEs:
\begin{align*} x_1'&=v_{1x}\\ v_{1x}'&=\frac{gm_2(x_2-x_1)}{((x_2-x_1)^2+(y_2-y_1)^2+(z_2-z_1)^2)^{\frac{3}{2}}}+\frac{gm_3(x_3-x_1)}{((x_3-x_1)^2+(y_3-y_1)^2+(z_3-z_1)^2)^{\frac{3}{2}}}\\ y_1'&=v_{1y}\\ v_{1y}'&=\frac{gm_2(y_2-y_1)}{((x_2-x_1)^2+(y_2-y_1)^2+(z_2-z_1)^2)^{\frac{3}{2}}}+\frac{gm_3(y_3-y_1)}{((x_3-x_1)^2+(y_3-y_1)^2+(z_3-z_1)^2)^{\frac{3}{2}}}\\ z_1'&=v_{1z}\\ v_{1z}'&=\frac{gm_2(z_2-z_1)}{((x_2-x_1)^2+(y_2-y_1)^2+(z_2-z_1)^2)^{\frac{3}{2}}}+\frac{gm_3(z_3-z_1)}{((x_3-x_1)^2+(y_3-y_1)^2+(z_3-z_1)^2)^{\frac{3}{2}}}\\ x_2'&=v_{2x}\\ v_{2x}'&=\frac{gm_1(x_1-x_2)}{((x_2-x_1)^2+(y_2-y_1)^2+(z_2-z_1)^2)^{\frac{3}{2}}}+\frac{gm_3(x_3-x_2)}{((x_3-x_2)^2+(y_3-y_2)^2+(z_3-z_2)^2)^{\frac{3}{2}}}\\ y_2'&=v_{2y}\\ v_{2y}'&=\frac{gm_1(y_1-y_2)}{((x_2-x_1)^2+(y_2-y_1)^2+(z_2-z_1)^2)^{\frac{3}{2}}}+\frac{gm_3(y_3-y_2)}{((x_3-x_2)^2+(y_3-y_2)^2+(z_3-z_2)^2)^{\frac{3}{2}}}\\ z_2'&=v_{2z}\\ v_{2z}'&=\frac{gm_1(z_1-z_2)}{((x_2-x_1)^2+(y_2-y_1)^2+(z_2-z_1)^2)^{\frac{3}{2}}}+\frac{gm_3(z_3-z_2)}{((x_3-x_2)^2+(y_3-y_2)^2+(z_3-z_2)^2)^{\frac{3}{2}}}\\ x_3'&=v_{3x}\\ v_{3x}'&=\frac{gm_1(x_1-x_3)}{((x_3-x_1)^2+(y_3-y_1)^2+(z_3-z_1)^2)^{\frac{3}{2}}}+\frac{gm_2(x_2-x_3)}{((x_3-x_2)^2+(y_3-y_2)+(z_3-z_2)^2)^{\frac{3}{2}}}\\ y_3'&=v_{3y}\\ v_{3y}'&=\frac{gm_1(y_1-y_3)}{((x_3-x_1)^2+(y_3-y_1)^2+(z_3-z_1)^2)^{\frac{3}{2}}}+\frac{gm_2(y_2-y_3)}{((x_3-x_2)^2+(y_3-y_2)^2+(z_3-z_2)^{\frac{3}{2}}}\\ z_3'&=v_{3z}\\ v_{3z}'&=\frac{gm_1(z_1-z_3)}{((x_3-x_1)^2+(y_3-y_1)^2+(z_3-z_1)^2)^{\frac{3}{2}}}+\frac{gm_2(z_2-z_3)}{((x_3-x_2)^2+(y_3-y_2)^2+(z_3-z_2)^{\frac{3}{2}}} \end{align*}
An example of a three-body problem in 3D, is setting the masses to \(m_1 = 0.03\), \(m_2 = 0.3\), and \(m_3 = 0.03\) and plotted the trajectories with these initial conditions:
\begin{align*} (x_1,y_1,z_1) &= (2,2,0) \\ (x_2,y_2,z_2) &= (0,0,0) \\ (x_3,y_3,z_3) &= (-2,0,-2) \\ (v_{x_1},v_{y_1},v_{z_1}) &= (0.2,-0.2,0) \\ (v_{x_2},v_{y_2},v_{z_2}) &= (0,0.01,-0.01) \\ (v_{x_3},v_{y_3},v_{z_3}) &= (-0.2,0.02,0.2) \end{align*}
Our code for the three-body problem in 3D is very similar to the part 2 and 3 code, just additional equations were added, as well as initial conditions. The video below shows the animation for this example.