  The strong magnetic field of neutron stars is intimately coupled to the
observed temperature and spectral properties, as well as to the observed timing
properties (distribution of spin periods and period derivatives). Thus, a
proper theoretical and numerical study of the magnetic field evolution
equations, supplemented with detailed calculations of microphysical properties
(heat and electrical conductivity, neutrino emission rates) is crucial to
understand how the strength and topology of the magnetic field vary as a
function of age, which in turn is the key to decipher the physical processes
behind the varied neutron star phenomenology. In this review, we go through the
basic theory describing the magneto-thermal evolution models of neutron stars,
focusing on numerical techniques, and providing a battery of benchmark tests to
be used as a reference for present and future code developments. We summarize
well-known results from axisymmetric cases, give a new look at the latest 3D
advances, and present an overview of the expectations for the field in the
coming years.
