The impact of strong magnetic fields B > 10¹³ G on the thermal evolution of neutron stars is investigated, including crustal heating by magnetic field decay. For this purpose, we perform 2D cooling simulations with anisotropic thermal conductivity considering all relevant neutrino emission processes for realistic neutron stars. The standard cooling models of neutron stars are called into question by showing that the magnetic field has relevant (and in many cases dominant) effects on the thermal evolution. The presence of the magnetic field significantly affects the thermal surface distribution and the cooling history of these objects during both the early neutrino cooling era and the late photon cooling era. The minimal cooling scenario is thus more complex than generally assumed. A consistent magnetothermal evolution of magnetized neutron stars is needed to explain the observations.