Calculations of energy balance in model stellar atmospheres for cold and moderately cold stars show that the overall energy gained or lost at the Lyman frequencies (lines plus continuum) can be up to 20 orders of magnitude higher than in the rest of the spectrum. Due to the great optical depth at those frequencies, however, monochromatic mean intensity of the radiation and the corresponding source function are numerically equal with many significant figures, and give therefore a vanishing net contribution to the energy balance. As all methods of temperature correction compute both integrals independently, this can lead to a biased temperature correction to the model output. We present here a method to calculate the correction of the temperature in stellar atmospheres that overcomes the aforesaid difficulty.