The topic of this dissertation is the fusion of a novel integration method, Modified Chebyshev Picard Iteration (MCPI), with Gauss' Variational Equations using a set of Modified Equinoctial Orbital Elements. This combination leads to a dramatically increased domain of Picard iteration convergence and an efficiency increase for MCPI solutions of the Initial Value Problem of Celestial Mechanics, thereby reducing the number of full gravity function calls. The set of Modified Equinoctial Orbital Elements (MEEs) are nonsingular over a large orbit variation domain, in contrast with the Classical Orbital Elements (COEs), which are singular at zero inclination and zero eccentricity, and propagation of MEEs with MCPI leads to much greater convergence time intervals for the IVP than is possible using Cartesian coordinates. This set of elements is also used to formulate the Two-Point Boundary Value Problem (TPBVP) associated with orbit transfer using a low-thrust, minimum-time control formulation and solves iteratively via a shooting method known as the Method of Particular Solutions.