A fundamental mathematical algorithm is presented for solving the wave equation in inhomogeneous media. This method completely generalizes the Haskell matrix method, which is the standard method for solving the wave equation in laterally homogeneous media. The Haskell matrix method has been the mathematical basis for many seismic techniques in exploration geophysics.In the method presented the medium is divided into layers and vertically averaged within each layer. The wave equation, within a layer, is then decoupled into an eigenvalue equation of the horizontal coordinate and a wave equation of the vertical coordinate. The eigenvalue equation is solved numerically. The vertical equation is solved analytically, once the eigenvalues are found. The solution throughout the medium is constructed by matching layer solutions at layer interfaces.The solution process of this method is 'modular,' in the sense that each layer corresponds to an independent module and all the modules together form the final, total solution. Such a modular solution process has the following advantages. First, in a 2-D problem, for example, each module is a 1-D problem, which is a much simpler problem numerically than the original full 2-D problem. Second, the module solutions can be used repeatedly to form the solution corresponding to different problems. For example, in modeling, only those layers which differ between two models require recalculation.The solution to plane-wave diffraction by a cylinder is obtained using this method, and it agrees well with the analytical solution.