Magneto-thermal evolutionary models are fundamental to explaining the phenomenological diversity of the different classes of neutron stars (NS) that compose the NS zoo. Although in the past decade the 2D magneto-thermal codes have been able to describe the main features of the cooling of highly magnetized NSs (see Pons & Viganò 2019 for a review), we know that to account for the non-axisymmetric modes arising from the magnetic field evolution (even for axisymmetric initial conditions) a 3D model is needed. Moreover, the thermal evolution in presence of a non-axisymmetric magnetic field leads to the formation of inhomogeneities in the NS surface temperature (e.g. hot spots), which can in principle explain the pulsed thermal emission often detected in highly magnetized NSs. Here I present the preliminary results of the thermal evolution part of the finite volume code MATINS (presented in Dehman et al. 2022 for its magnetic evolution part).