In this paper, we investigate the cooling of neutron stars with relativistic and nonrelativistic models of dense nuclear matter. We focus on the effects of uncertainties originated from the nuclear models, the composition of elements in the envelope region, and the formation of superfluidity in the core and the crust of neutron stars. Discovery of [Formula: see text] neutron stars PSR J1614−2230 and PSR J0343[Formula: see text]0432 has triggered the revival of stiff nuclear equation of state at high densities. In the meantime, observation of a neutron star in Cassiopeia A for more than 10 years has provided us with very accurate data for the thermal evolution of neutron stars. Both mass and temperature of neutron stars depend critically on the equation of state of nuclear matter, so we first search for nuclear models that satisfy the constraints from mass and temperature simultaneously within a reasonable range. With selected models, we explore the effects of element composition in the envelope region, and the existence of superfluidity in the core and the crust of neutron stars. Due to uncertainty in the composition of particles in the envelope region, we obtain a range of cooling curves that can cover substantial region of observation data.