New large observational surveys such as Gaia are leading us into an era of data abundance, offering unprecedented opportunities to discover new physical laws through the power of machine learning. Here we present an end-to-end strategy for recovering a free-form analytical potential from a mere snapshot of stellar positions and velocities. First we show how auto-differentiation can be used to capture an agnostic map of the gravitational potential and its underlying dark matter distribution in the form of a neural network. However, in the context of physics, neural networks are both a plague and a blessing as they are extremely flexible for modeling physical systems but largely consist in non-interpretable black boxes. Therefore, in addition, we show how a complementary symbolic regression approach can be used to open up this neural network into a physically meaningful expression. We demonstrate our strategy by recovering the potential of a toy isochrone system.