Once the need for an iterative procedure in order to solve the problem of the formation of spectral lines in the case of a model atom with many energy levels, the sequel is to seek for the most effective form of such an iterative scheme. It is an almost universal is assumed within all the iterative methods for the solution of those radiative transfer problems, in which the transfer equations are coupled to the state of the matter, to take as the input of each step of iterations the values of the opacity coefficients obtained as a result of the previous one. This is, for instance, the procedure used to correct the temperature in the computation of stellar atmosphere models, or that to build the -operator (either the exact or the approximated one) within the Accelerated Lambda Iteration methods. Yet, if we assume, in order to solve the multilevel line transfer problem, that at each step of iterations the line opacities are known, we can express via the statistical equilibrium equations the populations of the energy levels - and consequently the source functions of the relevant spectral lines - as a linear function of the full set of the corresponding mean intensities of the radiation field. Once such linear forms for the source functions, we are able to solve without the need of any further approximation the radiative transfer equations for are obtained lines, now linearly coupled through the above linear forms of the statistical equilibrium equations. This is achieved by means of the Implicit Integral Method that we already presented in a series of previous papers.