Modeling by paraxial extrapolators is applicable to wave-propagation problems in which most of the energy is traveling within a restricted angular cone about a principal axis of the problem. Using this technique, frequency-domain finite-difference solutions accurate for propagation angles out to 60 degrees are readily generated for both two-dimensional (2-D) and three-dimensional (3-D) models. Solutions for 3-D problems are computed by applying the 2-D paraxial operators twice, once along the x-axis and once along the y-axis, at each extrapolation step. The azimuthal anisotropy inherent to this splitting technique is essentially eliminated by adding a phase-correction operator to the extrapolation system. For heterogeneous models, scattering effects are incorporated by determining transmission and reflection coefficients at structural boundaries within the media. The direct forward-scattered waves are modeled with a single pass of the extrapolation operator in the paraxial direction for each frequency. The first-order backscattered energy is then modeled by extrapolation (in the opposite direction) of the reflected field determined on the first pass. Higher order scattering can be included by sweeping through the model with more passes.The chief advantages of the paraxial approach are (1) active storage is reduced by one dimension compared to solutions which must track both forward-scattered and backscattered waves simultaneously; thus, realistic 3-D problems can fit on today's computers, (2) the decomposition in frequency allows the technique to be implemented on highly parallel machines, (3) attenuation can be modeled as an arbitrary function of frequency, and (4) only a small number of frequencies are needed to produce movie-like time slices.