We present an efficient and accurate method for simulating massive neutrinos in cosmological structure formation simulations, together with an easy to use public implementation. Our method builds on our earlier implementation of the linear response approximation (LRA) for neutrinos, coupled with an N-body code for cold dark matter particles. The LRA's good behaviour at early times and in the linear regime is preserved, while better following the non-linear clustering of neutrinos on small scales. Massive neutrinos are split into initially "fast" and "slow" components. The fast component is followed analytically with the LRA all the way to redshift zero. The slow component is evolved with the LRA only down to a switch-on redshift $z_\nu = 1$, below which it is followed with the particle method, in order to fully account for its non-linear evolution. The slow neutrino particles are initialized at $z = 99$ in order to have accurate positions and velocities at the switch-on time, but are not used to compute the potential until $z \leq 1$, thus avoiding the worst effect of particle shot noise. We show that our hybrid method matches (and for small neutrino masses, exceeds) the accuracy of neutrino particle simulations with substantially lower particle load requirements.