The development of fast numerical methods for multilevel radiative transfer (RT) applications often leads to important breakthroughs in astrophysics, because they allow the investigation of problems that could not be properly tackled using the methods previously available. Probably, the most familiar example is the so-called Multilevel Accelerated Λ-Iteration (MALI) technique of Rybicki & Hummer for the case of a local approximate operator, which is based on Jacobi iteration. However, there are superior operator-splitting methods, based on Gauss-Seidel (GS) and Successive Overrelaxation (SOR) iteration, which provide a dramatic increase in the speed with which non-LTE multilevel transfer problems can be solved in one, two and three-dimensional geometries. Such RT methods, which were introduced by Trujillo Bueno & Fabiani Bendicho ten years ago, are the main subject of the first part of this paper. We show in some detail how they can be applied for solving multilevel RT problems in spherical geometry, for both atomic and molecular line transitions. The second part of the article addresses the issue of the calculation of the molecular number densities when the approximation of instantaneous chemical equilibrium turns out to be inadequate, which happens to be the case whenever the dynamical time scales of the astrophysical plasma under consideration are much shorter than the time needed by the molecules to form.