Boltzmann solvers are an important tool for the computation of cosmological observables in the linear regime. They involve solving the Boltzmann equation, followed by an integration in momentum space, to arrive at the desired fluid properties. This is a cumbersome, computationally expensive procedure. In this work we introduce the so-called generalized Boltzmann hierarchy (GBH) for massive neutrinos in cosmology, a simpler alternative to the usual Boltzmann hierarchy, where the momentum dependence is integrated out leaving us with a two-parameter infinite set of ordinary differential equations. Along with the usual expansion in multipoles, there is now also an expansion in higher velocity weight integrals of the distribution function. We show that the GBH produces the density contrast neutrino transfer function to a per mille level accuracy at both large and intermediate scales compared to the neutrino free-streaming scale. Furthermore, by introducing a switch to a viscous fluid approximation after horizon crossing, we show that the GBH can achieve over all scales the same accuracy as the standard CLASS approach in its default precision settings. The GBH is then a powerful tool to include neutrino anisotropies in the computation of cosmological observables in linear theory, with integration being simpler and potentially faster than standard methods.