An efficient recursive algorithm, the stiffness matrix method, has been developed for wave propagation in multilayered generally anisotropic media. This algorithm has the computational efficiency and simplicity of the standard transfer matrix method and is unconditionally computationally stable for high frequency and layer thickness. In this algorithm, the stiffness (compliance) matrix is calculated for each layer and recursively applied to generate a stiffness (compliance) matrix for a layered system. Next, reflection and transmission coefficients are calculated for layered media bounded by liquid or solid semispaces. The results show that the method is stable for arbitrary number and thickness of layers and the computation time is proportional to the number of layers. It is shown both numerically and analytically that for a thick structure the solution approaches the solution for a semispace. This algorithm is easily adaptable to laminates with periodicity, such as multiangle lay-up composites. The repetition and symmetry of the unit cell are naturally incorporated in the recursive scheme. As an example the angle beam time domain pulse reflections from fluid-loaded multilayered composites have been computed and compared with experiment. Based on this method, characteristic equations for Lamb waves and Floquet waves in periodic media have also been determined.