Maxwell's equations for propagation of electromagnetic waves in dispersive and absorptive (passive) media are represented in the form of the Schrödinger equation i∂Ψ/∂t = HΨ, where H is a linear differential operator (Hamiltonian) acting on a multi-dimensional vector Ψ composed of the electromagnetic fields and auxiliary matter fields describing the medium response. In this representation, the initial value problem is solved by applying the fundamental solution exp(−itH) to the initial field configuration. The Faber polynomial approximation of the fundamental solution is used to develop a numerical algorithm for propagation of broad band wave packets in passive media. The action of the Hamiltonian on the wave function Ψ is approximated by the Fourier grid pseudospectral method. The algorithm is global in time, meaning that the entire propagation can be carried out in just a few time steps. A typical time step is much larger than that in finite differencing schemes, ∆t F ≫ H −1 . The accuracy and stability of the algorithm is analyzed. The Faber propagation method is compared with the Lanczos-Arnoldi propagation method with an example of scattering of broad band laser pulses on a periodic grating made of a dielectric whose dispersive properties are described by the Rocard-Powels-Debye model. The Faber algorithm is shown to be more efficient. The Courant limit for time stepping, ∆t C ∼ H −1 , is exceeded at least in 3000 times in the Faber propagation scheme.