How to estimate which astronaut ends up furthest from the ISS after one orbit?

Question

Three astronauts; Goofus, Gallant and Zippy the Pinhead play a game. All three take their cold gas thruster-powered jet packs and quickly move 100 meters away from the ISS. For some reason after accelerating and decelerating to a stop, all three ran out of propellant.

Goofus moved 100 meters out in front of the ISS while Gallant moved to the side; and Zippy the Pinhead moved 100 meters farther from Earth so that he could be seen from ISS' purported zenith-facing window (1, 2, 3).

The four position vectors at the beginning for ISS, Goofy, Gallant and Zippy are:

x0 = [a, 0, 0.] + [a, 100, 0.] + [a, 0, 100.] + [a+100, 0, 0.]

and all four starting velocity vectors are [0, speed, 0.]

After one orbit around the Earth had elapsed, how far did each of our three heroes end up from the ISS?

---

This is what I got, but I'd prefer to accept an answer based more on math and principles than brute-force numerical integration.

J2 is important but for the purposes of other answers it can be ignored.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint

def deriv(X, t):
    x, v = X.reshape(2, -1)
    accs = []
    for xx in x.reshape(-1, 3):
        acc = -xx * GMe * ((xx**2).sum())**-1.5
        if use_J2:
            x, y, z = xx
            xsq, ysq, zsq = xx**2
            rm7 = (xsq + ysq + zsq)**-3.5
            accJ2x = x * rm7 * (6*zsq - 1.5*(xsq + ysq))
            accJ2y = y * rm7 * (6*zsq - 1.5*(xsq + ysq))
            accJ2z = z * rm7 * (3*zsq - 4.5*(xsq + ysq))
            acc -= J2 * np.hstack((accJ2x, accJ2y, accJ2z))
        accs.append(acc)
    return np.hstack((v, np.hstack((accs))))

halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]
GMe = 3.98600435436E+14 # earth (from DE430)
Re = 6378136.3 # meters
J2_earth = -1.08262545E-03  # unitless
J2 = J2_earth * (GMe * Re**2) 

inc = 51*pi/180.

a = Re + 400*1000.
speed = np.sqrt(GMe/a) # close enough for government work
x0 = [a, 0, 0.] + [a, 100, 0.] + [a, 0, 100.] + [a+100, 0, 0.]
v0 = 4*[0, speed, 0.]
X0 = np.array(x0 + v0, dtype=float).reshape(-1, 3)
s, c = [f(inc) for f in (np.sin, np.cos)]
R = np.array([[1, 0, 0], [0, c, -s], [0, s, c]])
X0 = (X0[:, None, :] * R).sum(axis=2).flatten()

times = np.arange(0, 92*60, 10)

use_J2 = True
answer, info = ODEint(deriv, X0, times, full_output=True)
use_J2 = False
answer_noJ2, info = ODEint(deriv, X0, times, full_output=True)

s0, s1 = answer.shape
answers = answer.T[:s1>>1].reshape(-1, 3, s0)
if True:
    plt.figure()
    for a in answers[1:]:
        r = np.sqrt(((a-answers[0])**2).sum(axis=0))
        plt.plot(times/60., r)
    plt.title('Goofus, Gallant and Zippy the Pinhead start 100 meters from ISS')
    plt.ylim(0, None)
    plt.xlabel('minutes')
    plt.show()