Spectral synthesis calculations in stellar (magnetized) atmospheres are based on the solution of the radiative transfer equation (RTE) for polarized light. The thermodynamic and magnetic properties of the atmospheres, along with the radiation field, completely specify the basic ingredients of the RTE, after which numerical methods have to be employed to calculate the emergent Stokes spectra. The advent of powerful analysis techniques for the inversion of Stokes spectra has evidenced the need for accurate and fast solutions of the RTE. In this paper we describe a novel Hermitian strategy to integrate the polarized RTE that is based on the Taylor expansion of the Stokes parameter vector to fourth order in depth. Our technique makes use of the first derivatives of the absorption matrix and source vector with respect to the coordinate measured along the ray path. Both analytical and numerical results indicate that the new strategy is superior to other methods in terms of speed and accuracy. It also gives an approximation to the evolution operator at no extra cost, which is of interest for inversion algorithms based on response functions. The Hermitian technique can be straightforwardly particularized to the scalar case, providing a very efficient solution of the RTE in the absence of magnetic fields. We investigate in detail the consequences of the oscillations that appear in the evolution operator for large values of line strength eta_0. The problems they pose are shared by all integration schemes, but can be minimized by adopting nonequally spaced grids.