  The thermal evolution of isolated neutron stars is a key element in
unraveling their internal structure and composition and establishing
evolutionary connections among different observational subclasses. Previous
studies have predominantly focused on one-dimensional or axisymmetric
two-dimensional models. In this study, we present the thermal evolution
component of the novel three-dimensional magnetothermal code MATINS
(MAgneto-Thermal evolution of Isolated Neutron Star). MATINS employs a finite
volume scheme and integrates a realistic background structure, along with
state-of-the-art microphysical calculations for the conductivities, neutrino
emissivities, heat capacity, and superfluid gap models. This paper outlines the
methodology employed to solve the thermal evolution equations in MATINS, along
with the microphysical implementation which is essential for the thermal
component. We test the accuracy of the code and present simulations with
non-evolving magnetic fields of different topologies to produce temperature
maps of the neutron star surface. Additionally, for a specific magnetic field
topology, we show one fully coupled evolution of magnetic field and
temperature. Subsequently, we use a ray-tracing code to link the neutron star
surface temperature maps obtained by MATINS with the phase-resolved spectra and
pulsed profiles that would be detected by distant observers. This study,
together with our previous article focused on the magnetic formalism, presents
in detail the most advanced evolutionary code for isolated neutron stars, with
the aim of comparison with their timing properties, thermal luminosities and
the associated X-ray light curves.
