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.