  We present in this work a novel and yet extremely simple method for
incorporating the effects of massive neutrinos in cosmological $N$-body
simulations. This so-called "SuperEasy linear response" approach is based upon
analytical solutions to the collisionless Boltzmann equation in the clustering
and free-streaming limits, which are then connected by a rational function
interpolation function with cosmology-dependent coefficients given by simple
algebraic expressions of the cosmological model parameters. The outcome is a
{\it one-line modification} to the gravitational potential that requires only
the cold matter density contrast as a real-time input, and that can be
incorporated into any $N$-body code with a Particle--Mesh component with no
additional implementation cost. To demonstrate its power, we implement the
SuperEasy method in the publicly available \gadgetcode{} code, and show that
for neutrino mass sums not exceeding $\sum m_\nu \simeq 1$ eV, the total matter
and cold matter power spectra are in sub-1\% and sub-0.1\% agreement with those
from state-of-the-art linear response simulations in literature. Aside from its
minimal implementation cost, compared with existing massive neutrino simulation
methods, the SuperEasy approach has better memory efficiency, incurs no runtime
overhead relative to a standard $\Lambda$CDM simulation, and requires no
post-processing. The minimal nature of the method allows limited computational
resources to be diverted to modelling other physical effects of interest, e.g.,
baryonic physics via hydrodynamics.
