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 . (read more)