Natural and engineered media usually involve combinations of solid, fluid and porous layers, and accurate and stable modelling of wave propagation in such complex multilayered media is fundamental to evaluating their properties with wave-based methods. Here we present a general stiffness matrix method for modelling waves in arbitrary multilayers. The method first formulates stiffness matrices for individual layers based on the governing wave equations for fluids and solids, and the Biot theory for porous materials. Then it utilises the boundary conditions considered at layer interfaces to assemble the layer matrices into a global system of equations, to obtain solutions for reflection and transmission coefficients at any incidence. Its advantage over existing methods is manifested by its unconditional computational stability, and its validity is proved by experimental validations on single solid sheets, porous layers, and porous-solid-porous battery electrodes. This establishes a powerful theoretical platform that allows us to develop advanced wave-based methods to quantitatively characterise properties of the layers, especially for layers of porous materials.