(* Reconstructed successful vector-equation Wolfram Language handle for 007.
   Uses the exact published day-two initial states and adaptive NDSolveValue. *)
ClearAll[u, w, t, acceleration, vars0, vars1, q0, q1, v0, equations0,
  equations1, solution0, solution1, marks, separations, recorded, relativeErrors];

acceleration[z_] := Flatten@Table[
  Sum[
    If[i == j, {0, 0, 0},
      With[{delta = z[[3 j - 2 ;; 3 j]] - z[[3 i - 2 ;; 3 i]]},
        delta/(delta.delta)^(3/2)
      ]
    ],
    {j, 3}
  ],
  {i, 3}
];

vars0 = Array[u, 9];
vars1 = Array[w, 9];
q0 = Flatten[{{0., 0., -1.}, {0., 0., 0.}, {0., 0., 1.}}];
q1 = Flatten[{{1.*^-7, 0., -1.}, {0., 0., 0.}, {-1.*^-7, 0., 1.}}];
v0 = Flatten[{{0., -1., 0.}, {1., 1., 0.}, {-1., 0., 0.}}];

equations0 = Join[
  Table[vars0[[i]]''[t] == acceleration[Through[vars0[t]]][[i]], {i, 9}],
  Table[vars0[[i]][0] == q0[[i]], {i, 9}],
  Table[vars0[[i]]'[0] == v0[[i]], {i, 9}]
];
equations1 = Join[
  Table[vars1[[i]]''[t] == acceleration[Through[vars1[t]]][[i]], {i, 9}],
  Table[vars1[[i]][0] == q1[[i]], {i, 9}],
  Table[vars1[[i]]'[0] == v0[[i]], {i, 9}]
];

solution0 = NDSolveValue[
  equations0, vars0, {t, 0, 80},
  AccuracyGoal -> 11, PrecisionGoal -> 11, MaxSteps -> Infinity
];
solution1 = NDSolveValue[
  equations1, vars1, {t, 0, 80},
  AccuracyGoal -> 11, PrecisionGoal -> 11, MaxSteps -> Infinity
];

marks = {20, 40, 60, 80};
separations = Table[
  Max[Norm /@ Partition[Through[solution0[mark]] - Through[solution1[mark]], 3]],
  {mark, marks}
];
recorded = {
  4.052889459131838*^-6, 1.8619198312937714*^-4,
  6.051802989520643*^-3, 1.366248921603192*^-1
};
relativeErrors = Abs[separations - recorded]/recorded;

Print[AssociationThread[marks, separations]];
Print["relative errors: ", relativeErrors];
Print["claim holds: ", TrueQ[Max[relativeErrors] < 0.005]];