The onset of thermal natural convection in a horizontal multilayer system consisting of a homogeneous porous layer sandwiched between two fluid layers has been simulated by an accurate numerical method. The porous and fluid layers include uniform heat sources. Flow in the porous medium has been governed by Darcy–Brinkman’s law. On the other hand, the Navier–Stokes equations with the Boussinesq approximation have ruled over the clear fluid layers. The lower and upper rigid surfaces are assumed to be fixed at the equal temperatures T_L and T_U. The eigenvalues and eigenfunctions of the linear stability analysis have been solved by utilizing the compound matrix method (CMM). The CMM reaches accurate results in a very efficient manner. Moreover, the method removes the stiffness from the equations of the stability system. The results indicate that the onset of convection and the nature of convection cells depend on the relative depths of layers. It has been observed that the thickness of the lower fluid layer increases the critical Rayleigh number of the upper fluid layer and stabilizes it.