This paper introduces an efficient computational procedure for analyzing the propagation of harmonic waves in layered elastic media. This offers several advantages, including the ability to handle arbitrary frequencies, depths, and the number of layers above an elastic half-space, and efforts to follow dispersion curves and flag up possible singularities are investigated. While there are inherent limitations in terms of computational accuracy and capacity, this methodology is straightforward to implement for studying free or forced vibrations and obtaining relevant response data. We present computations of wavenumber dispersion diagrams, phase velocity plots, and response data in both the frequency and time domains. These computational results are provided for two example cases: plane strain and axisymmetry. Our methodology is grounded in a well-conditioned dynamic stiffness approach specifically tailored for deep-layered strata analysis. We introduce an innovative method for efficiently computing wavenumber dispersion curves. By tracking the slope of these curves, users can effectively manage continuation parameters. We illustrate this technique through numerical evidence of a layer resonance in a real-life case study characterized by a fold in the dispersion curves. Furthermore, this framework is particularly advantageous for engineers addressing problems related to ground-borne vibrations. It enables the analysis of phenomena such as zero group velocity (ZGV), where a singularity occurs, both in the frequency and time domains, shedding light on the unique characteristics of such cases. Given the reduced dimension of the problem, this formulation can considerably aid geophysicists and engineers in areas such as MASW or SASW techniques.
