Two polynomial expansions of the time-evolution superoperator to directly integrate Markovian Liouville-von Neumann (LvN) equations for quantum open systems, namely the Newton interpolation and the Faber approximation, are presented and critically compared. Details on the numerical implementation including error control, and on the performance of either method are given. In a first physical application, a damped harmonic oscillator is considered. Then, the Faber approximation is applied to compute a condensed phase absorption spectrum, for which a semianalytical expression is derived. Finally, even more general applications are discussed. In all applications considered here it is found that both the Newton and Faber integrators are fast, general, stable, and accurate.