  Precision measurements of the number of effective relativistic neutrino
species and the primordial element abundances require accurate theoretical
predictions for early Universe observables in the Standard Model and beyond.
Given the complexity of accurately modelling the thermal history of the early
Universe, in this work, we extend a previous method presented by the author to
obtain simple, fast and accurate early Universe thermodynamics. The method is
based upon the approximation that all relevant species can be described by
thermal equilibrium distribution functions characterized by a temperature and a
chemical potential. We apply the method to neutrino decoupling in the Standard
Model and find $N_{\rm eff}^{\rm SM} = 3.045$ -- a result in excellent
agreement with previous state-of-the-art calculations. We apply the method to
study the thermal history of the Universe in the presence of a very light
($1\,\text{eV}<m_\phi < 1\,\text{MeV}$) and weakly coupled ($\lambda \lesssim
10^{-9}$) neutrinophilic scalar. We find our results to be in excellent
agreement with the solution to the exact Liouville equation. Finally, we
release a code: NUDEC_BSM (available in both Mathematica and Python formats),
with which neutrino decoupling can be accurately and efficiently solved in the
Standard Model and beyond: https://github.com/MiguelEA/nudec_BSM .
