Several methods have been proposed for correcting the temperature distribution in the calculation of model stellar atmospheres. The usually efficient procedure of Avrett and Krook (1963) has been widely used, although it may be superseded in some cases by the Auer and Mihalas (1968) method. Recently, still another procedure was proposed (Simonneau and Crivellari, 1988) which deserves further consideration. For our analysis we selected a method due to Napier and Dodd (1974) because it admits a rather simple numerical implementation. Also, as we were primarily interested in such things as how fast the method leads to final values, whether these values are the correct ones, and the dependence of the results on the initial temperature distribution, we simply considered a RE and LTE atmosphere, only affected by H and H⁻ continuous absorption. A special procedure was devised for calculating the flux as accurately as necessary. Three models were calculated, from Tₑ=10000K and g=10000cm/seg², to solar conditions, and standard optical depth ranging from 1E-4 to 100 (standard optical depth was selected at λ2000, or λ5000, according to model). The temperature distribution was represented using a Chebyshev polinomial expansion, and using sixth order approximation, a flux constancy better than 3% was obtained; the inclusion of more than 6 polinomials does not assure a better result. The convergence is rather fast: only 4 or 5 iterations are necessary.