Summary: A fundamental solution is derived for time harmonic elastic waves originating from a point source and propagating in a restricted class of three-dimensional, unbounded heterogeneous media which have a Poisson ratio of \(0.25\) and elastic moduli that vary quadratically with respect to the depth co-ordinate. The first step in the solution procedure is to transform the displacement vector in the equations of dynamic equilibrium through scaling by the square root of the position-dependent shear modulus. The constraints generated through this procedure are satisfied by quadratic (in the depth co-ordinate) profiles of the elastic moduli. During the next step, a double Fourier transform with respect to the horizontal co-ordinates is applied to the dynamic equilibrium equations, which assume a form amenable to solution by a first order matrix differential equation system. This latter system is solved using a series expansion due to the presence of non-constant matrix coefficients. The last step in recovering the fundamental solution is inversion of the double Fourier transform. This is accomplished numerically through use of the FFT, because complexity of the first order system approach precludes analytic inversion. Finally, some numerical examples serve to illustrate the present methodology.