Title: Enabling Equation Learning with the Bayesian Model Evidence via systematic $R^2$-elimination. Abstract: Deep learning is a powerful method for tasks like predictions and classification but lacks interpretability and analytic access. Instead of fitting up to millions of parameters, an intriguing alternative for a wide range of problems would be to learn the governing equations from data. The resulting models would be concise, parameters could be interpreted, the model could adjust to shifts in data, and analytic analysis would allow for extra insights. Common challenges are model complexity identification, stable feature selection, expressivity, computational feasibility, and scarce data. In our work, the mentioned challenges are addressed by combining existing methods in a novel way. We choose multiple regression as a framework and argue that a large space of model equations is captured. For feature selection, we exploit the computationally cheap coefficient of determination ($R^2$) for a model elimination process in a semi-comprehensive search. Final model selection is achieved by exact values of the Bayesian model evidence with empirical priors, which is known to identify suitable model complexity without relying on mass data. Random polynomials, an epidemiological model, and the Lorenz system are used as examples. For the Lorenz system, which is particularly challenging due to its chaotic nature, we demonstrate the favourable performance of our approach to existing state-of-the-art like SINDy. 