Observations of magnetars and some of the high magnetic field pulsars have shown that their thermal luminosity is systematically higher than that of classical radio-pulsars, thus confirming the idea that magnetic fields are involved in their X-ray emission. Here we present the results of 2D simulations of the fully coupled evolution of temperature and magnetic field in neutron stars, including the state-of-the-art kinetic coefficients and, for the first time, the important effect of the Hall term. After gathering and thoroughly re-analysing in a consistent way all the best available data on isolated, thermally emitting neutron stars, we compare our theoretical models to a data sample of 40 sources. We find that our evolutionary models can explain the phenomenological diversity of magnetars, high-B radio-pulsars, and isolated nearby neutron stars by only varying their initial magnetic field, mass and envelope composition. Nearly all sources appear to follow the expectations of the standard theoretical models. Finally, we discuss the expected outburst rates and the evolutionary links between different classes. Our results constitute a major step towards the grand unification of the isolated neutron star zoo.