In electromagnetic wave propagation problems, it is usually necessary to calculate the field quantities over a wide band of frequencies. In this paper, we develop a computationally-efficient scheme, which combines the finite element method (FEM) with the Pade approximation procedure, to derive the power series expansion of the unknown solution vector in terms of the frequency. Explicit power series expressions of the matrix operator are obtained for boundary value problems that are defined, not only over bounded spatial domains, but also over unbounded domains truncated either by an absorbing boundary condition (ABC) or by a perfectly matched layer (PML). It is shown that the FEM matrix is always a polynomial function of the frequency variable, even with the ABC or PML mesh truncations. The coefficients of the power series expansion are obtained iteratively, and certain a priori estimates are derived for the radius of convergence of this series expansion. Finally, Pade approximants are utilized to extend the region of convergence of the power series, enabling us to cover the frequency band with a minimum number of LU decompositions. (C) 1999 Elsevier Science S.A. All rights reserved.