We present 2D simulations of the cooling of neutron stars with strong magnetic fields (B > 10¹³ 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.