  We study the mutual influence of thermal and magnetic evolution in a neutron
star's crust in axial symmetry. Taking into account realistic microphysical
inputs, we find the heat released by Joule effect consistent with the
circulation of currents in the crust, and we incorporate its effects in 2D
cooling calculations. We solve the induction equation numerically using a
hybrid method (spectral in angles, but a finite--differences scheme in the
radial direction), coupled to the thermal diffusion equation. We present the
first long term 2D simulations of the coupled magneto-thermal evolution of
neutron stars. This substantially improves previous works in which a very crude
approximation in at least one of the parts (thermal or magnetic diffusion) has
been adopted. Our results show that the feedback between Joule heating and
magnetic diffusion is strong, resulting in a faster dissipation of the stronger
fields during the first million years of a NS's life. As a consequence, all
neutron stars born with fields larger than a critical value (about 5 10^13 G)
reach similar field strengths (approximately 2-3 10^{13} G) at late times.
Irrespectively of the initial magnetic field strength, after $10^6$ years the
temperature becomes so low that the magnetic diffusion timescale becomes longer
than the typical ages of radio--pulsars, thus resulting in apparently no
dissipation of the field in old NS. We also confirm the strong correlation
between the magnetic field and the surface temperature of relatively young NSs
discussed in preliminary works. The effective temperature of models with strong
internal toroidal components are systematically higher than those of models
with purely poloidal fields, due to the additional energy reservoir stored in
the toroidal field that is gradually released as the field dissipates.
