One of the major goals of the field of Milky Way dynamics is to recover the gravitational potential field. Mapping the potential would allow us to determine the spatial distribution of matter - both baryonic and dark - throughout the Galaxy. Taking advantage in recent developments in machine learning, we present a novel method for determining the gravitational field from a snapshot of the phase-space positions of stars, based only on minimal physical assumptions. Using a sample of observed stellar phase-space positions, we first fit a smooth and differentiable - yet highly flexible - model of the distribution function. Using the collisionless Boltzmann equation, we then find the gravitational potential - represented by a feed-forward neural network - that renders this distribution function stationary. This method is far more flexible than previous parametric methods, which fit narrow classes of analytic models to the data. This is a promising approach to uncovering the density structure of the Milky Way, using rich datasets of stellar kinematics that will soon become available.