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 making hypotheses on the solution itself, namely the specific intensity I (τ; n) of the radiation field. On the contrary, integral methods require making hypotheses on the source function S(τ). It seems of course more reasonable to make hypotheses on the latter because one can expect that the run of S(τ) with depth is smoother than that of I (τ; n). In previous works we assumed a piecewise parabolic approximation for the source function, which warrants the continuity of S(τ) and its first derivative at each depth point. Here we impose the continuity of the second derivative S''(τ). In other words, we adopt a cubic spline representation to the source function, which highly stabilizes the numerical processes.