There are serious discrepancies among some of the recent neutron star cooling calculations by various groups. We have been investigating the possible source of these discrepancies. In this paper, we report our findings. We also report the preliminary result of our most recent cooling calculations without assuming an isothermal stellar evolution code. In this work, we used the currently existing best energy transport theories, as well as general relativity, both in thermodynamics and hydrodynamics.