Magnetars are young neutron stars characterized by high X-ray quiescent luminosities, outbursts, and, in some cases, sporadic giant flares. They are believed to be powered by ultra-strong magnetic fields but the diversity of their observed behaviors is not understood and made even more puzzling by the discovery of "low-field" magnetars. By performing long-term two-dimensional simulations of the coupled magnetic, thermal, and rotational evolution of neutron stars, we try to establish evolutionary links between the (apparently different) observed phenomenology. With this approach we can also estimate the evolution of magnetic stresses in the crust, which allows us to establish when starquakes occur. From our results, we can give a qualitative description of the main stages in the evolution of a neutron star, from its youth to its old age, and connect them with observations.