We present a new flexible, fast and accurate way to implement massive neutrinos, warm dark matter and any other non-cold dark matter relics in Boltzmann codes. For whatever analytical or numerical form of the phase-space distribution function, the optimal sampling In momentum space compatible with a given level of accuracy is automatically found by comparing quadrature methods. The perturbation integration is made even faster by switching to an approximate viscous fluid description inside the Hubble radius, which differs front previous approximations discussed in the literature. When adding one massive neutrino to the minimal cosmological model, CLASS becomes just 1.5 times slower, instead of about 5 times in other codes (for fixed accuracy requirements). We illustrate the flexibility of our approach by considering a few examples of standard or non-standard neutrinos, as well as warm dark matter models.