The three-body problem (3BP) formulated by I. Newton has inspired many great mathematicians like L. Euler (1707-1783), J.L. Lagrange (1736-1813), K.G.J. Jacobi (1804-1851), W.R. Hamilton (1805-1865) and J.H. Poincaré (1854-1912) to develop mathematical studies, methods and theories in an attempt to solve this problem. Although these concepts have been crucial to the advancement of the natural sciences, a closed form solution has yet to be found. To aid in solving this problem, a restricted form of the 3BP has been formulated in which the mass of the two primaries greatly exceeds the mass of the third body. Known as the circular restricted three-body problem, this system gives rise to three collinear and two equilateral equilibria referred to as libration points. These locations are surrounded by various families of periodic libration point orbits. The corresponding exotic trajectories are highly non-Keplerian which offer desirable characteristics that have revolutionised space mission design. Crucial to this effort has been the theory on hyperbolic invariant manifolds. These topological structures asymptotically arrive at (stable) or depart from (unstable) a selected target orbit by exploiting the natural dynamics of the system. A connection which joins stable and unstable manifold trajectories with no discrepancy in state space, constitutes to a free (natural) transfer mechanism. In the case that two different equilibria are connected, this path is called a heteroclinic connection. The potential existence of these solutions are of significant scientific interest and are known to exist in the planar case. In this research, the current theory is extended by providing a comprehensive understanding of the phase space of the (spatial) vertical Lyapunov (V-L) orbits and their associated manifolds including the potential existence of heteroclinic connections in the Earth-Moon system.The transition from the planar to the spatial case is associated with a rapid increase in complexity. To reduce the order of the problem, the two hyperbolic manifolds which are to be connected are integrated until they intersect a plane called a Poincaré section. These stopping conditions form a N-1-dimensional subset of the phase space. A completely new effort is aimed at investigating the influence of the orientation of this section on the state vector discrepancy that arises from the connection of hyperbolic manifolds emanating from collinear libration points L1 and L2. The aim of this study is to reveal novel transfer trajectories outside the solution space confined by the traditional locations and orientations of the Poincaré section. The path towards creating these insights consist of three stages, all of which employ a variety of numerical methods. As these steps are sequential, the individual chapters of this thesis cannot be considered as stand-alone work. First of all, the target orbits are generated through the refinement (Differential Correction) and extension (Numerical Continuation) of analytic approximations based on Perturbation Theory. The families of horizontal Lyapunov (H-L), halo, axial, and V-L are analysed to provide an overview of the range of possible solutions. The periodicity of each orbit is verified numerically and the trajectory is validated by studying the so-called monodromy matrix. Secondly, the manifolds are constructed using the eigenspace of this fundamental matrix. Once more a high level of generality is achieved, this time by analysing the sets of manifolds associated with three different families of target orbits at three distinct energy levels. Each member of each manifold is verified through studying deviations in the Jacobi's constant, and validated by using symmetries. Lastly, the stopping conditions for integration of these hyperbolic trajectories are varied to analyse the sensitivity with respect to the orientation of the Poincaré section.This research has revealed that the most optimal connections for V-L orbits are found outside of the traditional stopping conditions. The ramifications of this behaviour are extensive and are observed to arise from the fact that these hyperbolic trajectories curve behind the Moon. Moreover, the orthographic projections display interesting differences in global stability. The manifolds associated with the V-L family retain their shape, whereas those corresponding to the H-L and halo families behave in a chaotic way over time. This enables an accurate approximation over extended integration periods. In addition, the V-L orbits have revealed the largest eigenvalue moduli across all families in both equilibria which indicate a minimum time to unwind from the target orbit. In conclusion, these insights provide a fresh perspective on the range of possible solutions and phenomena that can be further explored. Moreover, new techniques have been employed to ensure the mathematical validity of libration point orbits and their hyperbolic invariant manifolds. As these structures are inherently difficult to verify, this effort forms a valuable contribution to the TU Delft Astrodynamic Toolbox (C++) for future research.