The number of effective relativistic neutrino species represents a fundamental probe of the thermal history of the early Universe, and as such of the Standard Model of Particle Physics. Traditional approaches to the process of neutrino decoupling are either very technical and computationally expensive, or assume that neutrinos decouple instantaneously. In this work, we aim to fill the gap between these two approaches by modeling neutrino decoupling in terms of two simple coupled differential equations for the electromagnetic and neutrino sector temperatures, in which all the relevant interactions are taken into account and which allows for a straightforward implementation of BSM species. Upon including finite temperature QED corrections we reach an accuracy on N eff in the SM of 0.01. We illustrate the usefulness of this approach to neutrino decoupling by considering, in a model independent manner, the impact of MeV thermal dark matter on N eff . We show that Planck rules out electrophilic and neutrinophilic thermal dark matter particles of m < 3.0 MeV at 95% CL regardless of their spin, and of their annihilation being s-wave or p-wave. We point out that thermal dark matter particles with non-negligible interactions with both electrons and neutrinos are more elusive to CMB observations than purely electrophilic or neutrinophilic ones. In addition, assisted by the accuracy of our approach, we show that CMB Stage-IV experiments will generically test thermal dark matter particles with m 15 MeV. We make publicly available the codes developed for this study at https://github.com/MiguelEA/ nudec_BSM.