An exact stiffness matrix method is presented to evaluate the dynamic response of a multi‐layered poroelastic medium due to time‐harmonic loads and fluid sources applied in the interior of the layered medium. The system under consideration consists of N layers of different properties and thickness overlying a homogeneous half‐plane or a rigid base. Fourier integral transform is used with respect to the x‐co‐ordinate and the formulation is presented in the frequency domain. Fourier transforms of average displacements of the solid matrix and pore pressure at layer interfaces are considered as the basic unknowns. Exact stiffness (impedance) matrices describing the relationship between generalized displacement and force vectors of a layer of finite thickness and a half‐plane are derived explicitly in the Fourier‐frequency space by using rigorous analytical solutions for Biot's elastodynamic theory for porous media. The global stiffness matrix and the force vector of a layered system is assembled by considering the continuity of tractions and fluid flow at layer interfaces. The numerical solution of the global equation system for discrete values of Fourier transform parameter together with the application of numerical quadrature to evaluate inverse Fourier transform integrals yield the solutions for poroelastic fields. Numerical results for displacements and stresses of a few layered systems and vertical impedance of a rigid strip bonded to layered poroelastic media are presented. The advantages of the present method when compared to existing approximate stiffness methods and other methods based on the determination of layer arbitrary coefficients are discussed.