We study the mutual influence of thermal and magnetic evolution in a neutron stars 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 NSs 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.