Radiative transfer (RT) problems in which the source function includes a scattering-like integral are typical two-points boundary problems. Their solution via differential equations implies to make hypotheses on the solution itself, namely the specific intensity I(tau;n) of the radiation field. On the contrary, integral methods require to make hypotheses on the source function S(tau). It looks of course more reasonable to make hypotheses on the latter because one can expect that the run of S(tau) with depth be smoother than that of I(tau;n). In previous works we assumed a piece-wise parabolic approximation for the source function, which warrants the continuity of S(tau) and its first derivative at each depth point. Here we impose the continuity of the second derivative S"(tau). In other words, we adopt a cubic spline representation to the source function, which highly stabilize the numerical processes.

A mayor problem that arises in the computation of stellar atmosphere models is the self consistent determination of the temperature distribution via the constraint of energy conservation. The energy balance includes the gains due to the absorption of radiation and the losses due to emission. It is well known that within each one of the two above integrals the part corresponding to spectral ranges whose opacity X(nu) is huge can overcome by many orders of magnitude the part that corresponds to the remaining frequencies. On the other hand, at those frequencies where X(nu) is very large, the mean intensity J(nu) of the radiation field shall be equal, up to many significant digits, to the source function S(nu) and consequently to the Planck function B(nu,T). Then their net share to the energy balance shall be null, albeit separately their contributions to the gain and loss integrals are the most important numerically. Thus the spectral range whose physical contribution to the overall balance is null will dominate numerically both sides of the energy balance equation, and consequently the errors on the determination of J(nu) and S(nu) at these frequencies will falsify the balance. It is possible to circumvent the numerical problem brought about by the foregoing circumstances by solving the radiative transfer equation for the variable I(n,nu) - S(nu), instead of the customary intensity I(n,nu). We present here a novel iterative algorithm, based on iteration factors already employed by us with success, which makes it possible a fast correction of the temperature by computing directly the difference between the radiative losses and gains at each step of the iterations.