The next generation large scale structure observation experiments loom over the horizon with ever increasing precision. Theoretical predictions on non-linear clustering are continuously improved to be in a position to take advantage of the incoming data. Deep within the non-linear regime, N-body simulations remain to be the only reliable method for calculating matter clustering. The immense computational resource cost of such simulations limit most to only cold dark matter. Concerted efforts in recent decades have been made to allow the inclusion of additional components, such as massive neutrinos. In this talk, I will be presenting two novel techniques to incorporate massive neutrinos as linear response perturbations into an N-body simulation, capable of achieving precent level accuracy for the matter power spectra, with small incremental cost to the computational requirements. These will serve as important steppingstones to next generation high resolution simulations with a more complete picture of the cosmos.