We present a model for calculating the evolution of neutron stars in the thermal evolution approximation, which follows the detailed temperature dependence in the interior of the star assuming a fixed density structure. A new numerical algorithm for solving the thermal evolution equations is presented and comparisons are made with other algorithms. Other numerical details of the model are discussed. The neutrino emissivities, heat capacities, thermal conductivities, and superfluid treatment used in the models are described. A comparison is made to the full-evolution models of Nomoto and Tsuruta. Results of the thermal evolution model are presented for representative neutron star models of a gravitational mass M = 1.4 M. with and without strong surface magnetic fields, and also for stars with a central accelerated cooling mechanism. Isothermal cooling models are examined to determine the ages beyond which the cooling model applies. Our predicted luminosities are compared with soft X-ray observations made with the Einstein Observatory.