A linear stability analysis is performed to investigate the stability of a plane Poiseuille flow of a fluid-porous-fluid three-layer system in a channel with smooth rigid walls. It is widely recognized that a channel with an isotropic porous layer exhibits more stability than those bounded by two porous walls. We investigate the specific role of anisotropy in this context by assuming an anisotropic porous layer saturated with the same fluid. The porous layer is modeled using the volume-averaged Navier–Stokes equation, and fluid layers using the Navier–Stokes equation. The flows are coupled through a stress-jump condition at the porous-fluid interface, allowing for a continuous velocity profile across the entire flow domain. The Orr–Sommerfeld type boundary value problem is developed to provide a framework for performing linear stability analysis. The spectral collocation method is used to solve the eigenvalue problem for the amplitude of arbitrary wavenumber disturbances. The instability characteristics originating from two-dimensional infinitesimal perturbations of the most unstable modes are investigated intensively. Three types of instability are identified in long-wave, moderate-wave, and short-wave regimes. The instability modes of the short-wave regime are specified as the fluid modes, and the modes in the other two regimes are recognized as the porous modes. Due to the presence of different modes, the neutral stability boundaries display uni-modal, bi-modal, and tri-modal structure that depends on the mean permeability and anisotropy. It is revealed that either an increase in mean permeability, a decrease in the anisotropy parameter, or a decrease in the depth ratio decreases the critical Reynolds number of each unstable mode, destabilizing the flow. At low wall normal permeability of the porous medium with large depth ratios (such as in the ground structure of geothermal systems), the porous mode is the dominant mode, triggering long-wave instability. Moreover, an energy budget analysis decodes the instability mechanism due to physical parameters. An equation for the rate of change of the kinetic energy of the two-dimensional infinitesimal disturbances is derived. Through the Reynolds stress, the energy production term transfers energy from the base flow to the disturbance, increasing the kinetic energies in all the layers and enhancing the growth rate of the least stable mode. The study paves the way for a stability analysis of a hydro-dynamically and thermally fully developed laminar flow in a multi-layer fluid-porous-fluid system that would facilitate the prediction of parameter regime where the overall performance of the systems can be optimized depending on the relevant applications.