The wave field is decomposed into its frequency-wave-number components, and a compound-matrix ODE theory is formulated for multiregion fluid–solid media using multipoint boundary conditions. Forward propagation of compound entities is done to a matchingdepth, from which stabilized backward propagation to the receivers is performed. Wave-number integration as well as modal synthesis is covered. Source arrays and receivers may be arbitrarily located. An appropriate adjointproblem is defined and solved. The formula for modal excitation coefficients is generalized to cover leakymodes, which have an exponential increase with depth. Three methods for stable and automatic computation of modal depth functions, in solid as well as fluid regions, are proposed. In work by M. B. Porter and E. L. Reiss [J. Acoust. Soc. Am. 77, 1760–1767 (1985)] the mode shapes in the solid bottom were left aside, and the method by F. Schwab etal. [Bull. Seismol. Soc. Am. 74, 1555–1578 (1984)] involves experimentation with a cutoff depth for an artificial homogeneous half-space. For a medium composed of homogeneous layers, it is shown how efficient computations are obtained by writing the propagator matrices involved as products of sparse matrices, which are applied in sequence. Numerical stability can be achieved without artificially splitting thick homogeneous layers.