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. 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. We first train a normalizing flow on a sample of observed phase-space positions, obtaining a smooth, differentiable approximation of the phase-space 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.