Teleseismic waves propagating in the upper mantle are subject to considerable distortion due to the effects of laterally heterogeneous structure. The magnitude and scale of velocity contrasts representative of features such as subducted slabs may be such that wave diffraction becomes an important process. In this case forward modeling methods based on high-frequency asymptotic approximations to the wave equation will not accurately describe the wavefield. A method is introduced to model the propagation of teleseismic P waves in a laterally heterogeneous upper mantle that accounts for distortion of the initial portion of the wavefield including the effects of multipathing and frequency-dependent diffraction. The method is based on a parabolic approximation to the wave equation that is solved in the time domain on a finite-difference grid which mimics the expected pattern of energy flow in a reference velocity field. Numerical examples for a simple two-dimensional subducting slab model demonstrate the application of the method and illustrate the effects of multipathing and diffraction which dominate waveform distortion at high and low frequencies, respectively.