  Spherical Jeans modeling is widely used to estimate mass profiles of systems
from star clusters to galactic stellar haloes to clusters of galaxies. It
derives the cumulative mass profile, M(<r), from kinematics of tracers of the
potential under the assumptions of spherical symmetry and dynamical
equilibrium. We consider the application of Jeans modeling to mapping the dark
matter distribution in the outer reaches of the Milky Way using field halo
stars. We present a novel non-parametric routine for solving the spherical
Jeans equation by fitting B-splines to the velocity and density profiles of
halo stars. While most implementations assume parametric forms for these
profiles, B-splines provide non-parametric fitting curves with analytical
derivatives. Our routine recovers the mass profiles of equilibrium systems with
flattened haloes or a stellar disc and bulge excellently (<~ 10% error at most
radii). Tests with non-equilibrium, Milky Way-like galaxies from the Latte
suite of FIRE-2 simulations perform quite well (<~ 15% error for r <~ 100 kpc).
We also create observationally motivated datasets for the Latte suite by
imposing selection functions and errors on phase space coordinates
characteristic of Gaia and the DESI Milky Way Survey. The resulting imprecise
and incomplete data require us to introduce an MCMC-based subroutine to obtain
deconvolved density and velocity dispersion profiles from the tracer
population. With these observational effects taken into account, the accuracy
of the Jeans mass estimate remains at the level 20% or better.
