In some recent publications, formulations of a time-domain beam propagation method (TD-BPM) were developed to incorporate dispersive media with complex permittivity described by a single first-order pole or a modified version thereof. However, for many emerging optical applications, procedures are needed for general linear dispersive media with dispersions described by rational functions, such as a combination of multiple second-order poles. The authors present an efficient frequency-dependent full-band TD-BPM formulation for the analysis of general linear dispersive media. In this method, the complex time-domain convolution of the dispersive media is evaluated efficiently via a Z-transform. The formulation is validated by simulating dispersive media with complex permittivity described by single and multiple second-order poles at microwave and optical frequencies. It is shown that the proposed formulation can accurately simulate broadband electromagnetic responses using a much larger time-step size than those required by other conventional numerical techniques. Extension of this formulation to the treatment of any general material, such as dispersive non-linear material, is straightforward.