Ever since the launch of the Luna 1 rover in 1959, scientists have noticed a peculiar phenomenon about Lunar surface gravity: as satellites and probes orbited the Moon over certain craters and impact basins, they periodically veered off course and sometimes even plummeted toward the lunar surface. The Moon’s surface gravity has always puzzled scientists in that its surface gravity anomaly poses many questions about the mass below the lunar surface. We propose a method utilizing orthogonal polynomial approximation in higher dimension spaces to model the surface gravity anomaly of the Moon. The Finite Element Method (FEM) is a numerical method for solving differential equations which can be seen in engineering, mathematical modeling, and other computationally demanding applications. With this method, we wish to approximate solutions to Laplace’s classical differential equations of gravity. Using classical Chebyshev polynomials as basis functions, the least-squares approximation was utilized to review discrete samples of the approximation function. The first application in this project was to substitute the classical spherical harmonic series of approximation with a series of orthogonal polynomial approximations (i.e using the FEM approach). With an error tolerance set at 10−9 m s−2, we use this method to adapt the gravity model radially upwards from the lunar surface to demonstrate how a higher degree order of approximation is required on and near the lunar surface, while a lower degree polynomial order is required as radius increases. We show that a five order of magnitude computational speedup is possible when applying the method to radial adaptation. More intrinsically, our second application involved using the methodology as an effective tool in solving boundary value problems. Specifically, we implement this approach to solve for classical differential equations involved with high precision, long-term orbit propagation. This application provided four order of magnitude speed up in computational time alongside a 10−10 m s−2 error for a series of orbit propagation tests. Alongside the advancements in orthogonal approximation theory, FEM enables revolutionary speedups in orbit propagation without loss in accuracy.