Heteroclinic connections in the spatial circular restricted three-body problem play an important role in astrodynamics. This paper presents a fully numerical methodology for the computation of the two-parameter families of connections that exist between families of invariant tori in the three-body problem. The computation of connections is presented as a two-point boundary value problem in which the initial and final states belong to unstable and stable manifolds of two respective normally hyperbolic invariant manifolds. A flow map torus computation procedure is modified to compute particular trajectories within these spaces, and a robust two-parameter continuation scheme is paired with a shooting technique to compute the entire connection family. Results are provided in the Sun–Earth and Earth–Moon representations of the spatial three-body problem. The Sun–Earth connection family works to verify the methodology with previously developed semi-analytical connection computation techniques. The Earth–Moon example demonstrates the methodology in a setting where fully numerical techniques are required to study the connection family.