Stable and accurate surface‐wave calculations present a long‐standing numerical and analytical challenge to seismologists. In a layered medium, we can describe the surface‐wave behavior with harmonic time dependence in terms of an eigenproblem where the eigenvalues and eigenfunctions define the surface‐wave modes.

Numerous studies have explored diverse numerical approaches to solve this eigenproblem, but they tend to suffer from numerical difficulties that limit the complexity of the medium, frequency range of applicability, or accuracy of the solution. We propose an equivalent formulation that replaces the conventional stress‐displacement vector with an alternative one, in order to cast the eigenproblem in a standard form that is linear in the eigenvalues. We discretize the system and boundary conditions using a Chebyshev spectral collocation method, leading to a finite‐dimensional generalized matrix eigenvalue problem that can be solved directly. No iterations are required to satisfy the boundary conditions or to isolate the eigenmodes. Collocation methods allow for the solution of the eigenproblem for general depth‐dependent elastic properties, including continuous depth variations of the properties as well as material interfaces. We illustrate the use of our technique to calculate dispersion curves, theoretical waveform time series, and to estimate the change from retrograde to prograde particle motion at the surface for a complex 1D structure under Los Angeles.