  We investigate the cooling of neutron stars with relativistic and
non-relativistic 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 $2 M_\odot$ neutron stars PSR J1614-2230
and PSR J0343+0432 has triggered the revival of stiff nuclear equation of state
at high densities. In the mean time, 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 evenlope 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.
