Our best estimates of the age, contents, and geometry of the Universe come from comparing predictions of the Einstein-Boltzmann (E-B) equations with observations of galaxies and the afterglow of the Big Bang. Existing E-B solvers are not differentiable, and Bayesian parameter estimation of these differential equation models are thus restricted to employing gradient-free inference algorithms. This becomes intractable in the high-dimensional settings increasingly relevant for modern observations. Propagating derivatives through the numerical solution of these ordinary differential equations is tractable through automatic differentiation (AD). We are actively developing the first AD-enabled E-B solver, Bolt.jl, making use of the rich Julia ecosystem of AD tools. Beyond mitigating the cost of high-dimensional inference, Bolt.jl opens the door to testing new cosmological physics against data at the level of terms in the Einstein-Boltzmann equations, using neural ODEs and physics-informed neural networks (PINNs).