We consider the impact of neutrino self-interactions described by an effective four-fermion coupling on cosmological observations. Implementing the exact Boltzmann hierarchy for interacting neutrinos first derived in  into the Boltzmann solver CLASS, we perform a detailed numerical analysis of the effects of the interaction on the cosmic microwave background (CMB) anisotropies, and compare our results with known approximations in the literature. While we find good agreement between our exact approach and the relaxation time approximation used in some recent studies, the popular (c2eff,c2vis)-parameterisation fails to reproduce the correct scale dependence of the CMB temperature power spectrum. We then proceed to derive constraints on the effective coupling constant Geff using currently available cosmological data via an MCMC analysis. Interestingly, our results reveal a bimodal posterior distribution, where one mode represents the standard ΛCDM limit with Geff≲108GF, and the other a scenario in which neutrinos self-interact with an effective coupling constant Geff≃3×109GF.