The thermal evolution of isothermal neutron stars is studied with matter both in the hadronic phase as well as in the mixed phase of hadronic matter and strange quark matter. In our models, the dominant early-stage cooling process is neutrino emission via the direct Urca process. As a consequence, the cooling curves fall too fast compared to observations. However, when superfluidity is included, the cooling of the neutron stars is significantly slowed down. Furthermore, we find that the cooling curves are not very sensitive to the precise details of the mixing between the hadronic phase and the quark phase and also of the pairing that leads to superfluidity.