The spherical Jeans equation is widely used to estimate mass profiles of systems from star clusters to galactic stellar halos to clusters of galaxies. The cumulative mass profile, M(<r), is derived from kinematics of tracers of the potential under the assumptions of spherical symmetry and dynamical equilibrium. We consider the application of the Jeans equation to mapping the outer reaches of the Milky Way, specifically to determine the dark matter distribution from field halo stars. We present a new non-parametric routine for solving the spherical Jeans equation by fitting B-Splines to the 3-dimensional velocity and density distributions of halo stars obtained by the Gaia survey and spectroscopic surveys such as DESI. While most implementations of the Jeans method assume parametric forms for these profiles, B-Splines provide non-parametric fitting curves and make no sacrifice in the convenience or accuracy of their derivatives. Despite Jeans Modeling’s prevalence, there is little work quantifying the errors introduced when breaking the assumptions of spherical symmetry and dynamical equilibrium. We validate our routine on several progressively more complex and realistic mock datasets that break these assumptions. We find that our spherical Jeans routine recovers the mass profiles of even quite flattened halos and systems including a stellar disk and bulge very well (<5% error). However, our tests with a mock data set from the cosmological hydrodynamic Auriga simulations (galaxy 6), yield significantly larger errors on the mass profile (~30%). This larger error suggests that the output of spherical Jeans modeling is more sensitive to deviations from dynamical equilibrium and the presence of substructure in the halo than deviations from sphericity. This work is funded in part by NASA grants NNX15AK79G and 80NSSC20K0509 and a MICDE Catalyst Grant.