  Simulating the long-term evolution of temperature and magnetic fields in
neutron stars is a major effort in astrophysics, having significant impact in
several topics. A detailed evolutionary model requires, at the same time, the
numerical solution of the heat diffusion equation, the use of appropriate
numerical methods to control non-linear terms in the induction equation, and
the local calculation of realistic microphysics coefficients. Here we present
the latest extension of the magneto-thermal 2D code in which we have coupled
the crustal evolution to the core evolution, including ambipolar diffusion. It
has also gained in modularity, accuracy, and efficiency. We revise the most
suitable numerical methods to accurately simulate magnetar-like magnetic
fields, reproducing the Hall-driven magnetic discontinuities. From the point of
view of computational performance, most of the load falls on the calculation of
microphysics coefficients. To a lesser extent, the thermal evolution part is
also computationally expensive because it requires large matrix inversions due
to the use of an implicit method. We show two representative case studies: (i)
a non-trivial multipolar configuration confined to the crust, displaying
long-lived small-scale structures and discontinuities; and (ii) a preliminary
study of ambipolar diffusion in normal matter. The latter acts on timescales
that are too long to have relevant effects on the timescales of interest but
sets the stage for future works where superfluid and superconductivity need to
be included.
