  We present 2D simulations of the cooling of neutron stars with strong
magnetic fields (B \geq 10^{13} G). We solve the diffusion equation in axial
symmetry including the state of the art microphysics that controls the cooling
such as slow/fast neutrino processes, superfluidity, as well as possible
heating mechanisms. We study how the cooling curves depend on the the magnetic
field strength and geometry. Special attention is given to discuss the
influence of magnetic field decay. We show that Joule heating effects are very
large and in some cases control the thermal evolution. We characterize the
temperature anisotropy induced by the magnetic field for the early and late
stages of the evolution of isolated neutron stars.
