Three-Body Problem in 3D

G = 1;
time = 20;
spScale = 8;

mA = 1.3;
xA0 = 0;
yA0 = 0;
zA0 = 0;
vxA0 = 0;
vyA0 = 0;
vzA0 = -0.5;

mB = 0.9;
xB0 = 1;
yB0 = 0;
zB0 = 0;
vxB0 = 0;
vyB0 = 0.8;
vzB0 = -0;

mC = 0.1;
xC0 = -1;
yC0 = 0;
zC0 = 0;
vxC0 = 0;
vyC0 = -0.8;
vzC0 = 0.2;

soln1 = NDSolve[
  {
   mA xA''[t] == -((
      G mA mB (xA[t] - 
      xB[t]))/((xA[t] - xB[t])^2 + (yA[t] - yB[t])^2 + (zA[t] - 
      zB[t])^2)^(3/2)) - (
      G mA mC (xA[t] - 
      xC[t]))/((xA[t] - xC[t])^2 + (yA[t] - yC[t])^2 + (zA[t] - 
      zC[t])^2)^(3/2),
   mA yA''[t] == -((
      G mA mB (yA[t] - 
      yB[t]))/((xA[t] - xB[t])^2 + (yA[t] - yB[t])^2 + (zA[t] - 
      zB[t])^2)^(3/2)) - (
      G mA mC (yA[t] - 
      yC[t]))/((xA[t] - xC[t])^2 + (yA[t] - yC[t])^2 + (zA[t] - 
      zC[t])^2)^(3/2),
   mA zA''[t] == -((
      G mA mB (zA[t] - 
      zB[t]))/((xA[t] - xB[t])^2 + (yA[t] - yB[t])^2 + (zA[t] - 
      zB[t])^2)^(3/2)) - (
      G mA mC (zA[t] - 
     zC[t]))/((xA[t] - xC[t])^2 + (yA[t] - yC[t])^2 + (zA[t] - 
     zC[t])^2)^(3/2),
   mB xB''[t] == -((
      G mB mC (xB[t] - 
      xC[t]))/((xA[t] - xB[t])^2 + (yA[t] - yB[t])^2 + (zA[t] - 
      zB[t])^2)^(3/2)) - (
      G mB mA (xB[t] - 
      xA[t]))/((xA[t] - xB[t])^2 + (yA[t] - yB[t])^2 + (zA[t] - 
      zB[t])^2)^(3/2),
   mB yB''[t] == -((
      G mB mC (yB[t] - 
      yC[t]))/((xA[t] - xB[t])^2 + (yA[t] - yB[t])^2 + (zA[t] - 
      zB[t])^2)^(3/2)) - (
      G mB mA (yB[t] - 
      yA[t]))/((xA[t] - xB[t])^2 + (yA[t] - yB[t])^2 + (zA[t] - 
      zB[t])^2)^(3/2),
   mB zB''[t] == -((
      G mB mC (zA[t] - 
      zC[t]))/((xA[t] - xB[t])^2 + (yA[t] - yB[t])^2 + (zA[t] - 
      zB[t])^2)^(3/2)) - (
      G mB mA (zB[t] - 
      zA[t]))/((xA[t] - xC[t])^2 + (yA[t] - yC[t])^2 + (zA[t] - 
      zC[t])^2)^(3/2),
   mC xC''[t] == -((
      G mC mA (xC[t] - 
      xA[t]))/((xA[t] - xB[t])^2 + (yA[t] - yB[t])^2 + (zA[t] - 
      zB[t])^2)^(3/2)) - (
      G mC mB (xC[t] - 
      xB[t]))/((xA[t] - xB[t])^2 + (yA[t] - yB[t])^2 + (zA[t] - 
      zB[t])^2)^(3/2),
   mC yC''[t] == -((
      G mC mA (yC[t] - 
      yA[t]))/((xA[t] - xB[t])^2 + (yA[t] - yB[t])^2 + (zA[t] - 
      zB[t])^2)^(3/2)) - (
      G mC mB (yC[t] - 
      yB[t]))/((xA[t] - xB[t])^2 + (yA[t] - yB[t])^2 + (zA[t] - 
      zB[t])^2)^(3/2),
   mB zC''[t] == -((
      G mC mA (zC[t] - 
      zA[t]))/((xA[t] - xB[t])^2 + (yA[t] - yB[t])^2 + (zA[t] - 
      zB[t])^2)^(3/2)) - (
      G mC mB (zC[t] - 
      zB [t]))/((xA[t] - xC[t])^2 + (yA[t] - yC[t])^2 + (zA[t] - 
      zC[t])^2)^(3/2),
   xA[0] == xA0, yA[0] == yA0, zA[0] == zA0, xA'[0] == vxA0, 
   yA'[0] == vyA0, zA'[0] == vzA0,
   xB[0] == xB0, yB[0] == yB0, zB[0] == zB0, xB'[0] == vxB0, 
   yB'[0] == vyB0, zB'[0] == vzB0,
   xC[0] == xC0, yC[0] == yC0, zC[0] == zC0, xC'[0] == vxC0, 
   yC'[0] == vyC0, zC'[0] == vzC0
   },
  {xA, yA, zA, xB, yB, zB, xC, yC, zC},
  {t, 0, time},
  MaxSteps -> 300000
  ]

x1[t_] := Evaluate[xA[t] /. soln1[[1, 1]]]
y1[t_] := Evaluate[yA[t] /. soln1[[1, 2]]]
z1[t_] := Evaluate[zA[t] /. soln1[[1, 3]]]
x2[t_] := Evaluate[xB[t] /. soln1[[1, 4]]]
y2[t_] := Evaluate[yB[t] /. soln1[[1, 5]]]
z2[t_] := Evaluate[zB[t] /. soln1[[1, 6]]]
x3[t_] := Evaluate[xC[t] /. soln1[[1, 7]]]
y3[t_] := Evaluate[yC[t] /. soln1[[1, 8]]]
z3[t_] := Evaluate[zC[t] /. soln1[[1, 9]]]

Manipulate[
 Show[
  {ParametricPlot3D[
    {
     {x1[t], y1[t], z1[t]},
     {x2[t], y2[t], z2[t]},
     {x3[t], y3[t], z3[t]}
     },
    {t, 1.005, tmax},
    Axes -> False,
    Ticks -> True,
    Boxed -> False,
    PlotRange -> {{-1.25, 0.75}, {0, 3.5}, {-0.5, 4}},
    PlotStyle -> {Lighter[Red, 0.5], Lighter[Green, 0.5], 
      Lighter[Blue, 0.5]},
    FaceGrids -> {{0, 0, -1}, {0, 1, 0}, {-1, 0, 0}}
    ]
   },
  {ParametricPlot3D[
    {
     {x1[t], y1[t], z1[t]},
     {x2[t], y2[t], z2[t]},
     {x3[t], y3[t], z3[t]}
     },
    {t, tmax - 1, tmax},
    PlotStyle -> {{Thick, Red}, {Thick, Green}, {Thick, Blue}}
    ]
   },
  {Graphics3D[
    {Red, Sphere[{x1[tmax], y1[tmax], z1[tmax]}, mA/spScale],
     Green, Sphere[{x2[tmax], y2[tmax], z2[tmax]}, mB/spScale],
     Blue, Sphere[{x3[tmax], y3[tmax], z3[tmax]}, mC/spScale]}
    ]
   },
  ImageSize -> {500, 500}
  ],
 {tmax, 1, 8.8}
 ]