  Efficient numerical approximation of the polarized radiative transfer
equation is challenging because this system of ordinary differential equations
exhibits stiff behavior, which potentially results in numerical instability.
This negatively impacts the accuracy of formal solvers, and small step-sizes
are often necessary to retrieve physical solutions. This work presents
stability analyses of formal solvers for the radiative transfer equation of
polarized light, identifies instability issues, and suggests practical
remedies. In particular, the assumptions and the limitations of the stability
analysis of Runge-Kutta methods play a crucial role. On this basis, a suitable
and pragmatic formal solver is outlined and tested. An insightful comparison to
the scalar radiative transfer equation is also presented.
