The design and simulation of ultrashort pulse shaping systems require pulse propagation methods that take the combined effects of dispersion, diffraction, and system aberrations into account. In the conventional pulse propagation methods based on the spectrum of plane wave, usually large number of sampling points are needed for the correct Fourier transform operations due to the fast oscillating phase of the complex pulse field. In this work, I have developed an alternative pulse propagation method, based on the Gaussian pulsed beam decomposition, as an extension of the monochromatic Gaussian beam propagation method. Methods for the decomposition of an input pulse, with arbitrary spatial and temporal (spectral) profiles, into a set of elementary Gaussian pulsed beams are proposed. Algorithms for computing the spatio-temporal and spatio-spectral profiles of the propagated pulse as the phase correct superposition of individual Gaussian pulsed beams are developed. The proposed decomposition method allows the elementary Gaussian pulsed beams to have different parameters depending on the local spatial and spectral phase of the given input pulse which reduces the number of Gaussian pulsed beams required to decompose an input pulsed beam with a given accuracy. Furthermore, a new kind of beam called the truncated Gaussian beam, is introduced and combined with the conventional Gaussian beam decomposition method to enable decomposition of fields after hard apertures. The analytical propagation equation of the truncated Gaussian beam through a paraxial optical system is derived. Additionally, the application of the Gaussian beam decomposition method is extended to handle the propagation of vectorial fields. Several example calculations are presented to validate the proposed methods and show their application in propagating fields through optical systems which are rather complicated to model using the conventional methods.