<p>This study investigates the pressure distribution in porous media while accounting for the memory effects in fluid flow. A mathematical model is developed to describe fluid flow in both homogeneous and heterogeneous porous media under various boundary conditions. The presence of porous media can lead to anomalous diffusion, which preserves the history of the interactions between the fluid and medium. To capture this behavior, the Caputo-type fractional derivative is incorporated into Darcy’s law. The problem is analyzed numerically using a spectral method based on orthogonal polynomials. In this approach, the derivatives of orthogonal polynomials are represented in matrix form, which is known as the operational matrices. These matrices provide a systematic framework for handling polynomials, facilitating their application to fluid flow problems. By employing operational matrices, the original fractional differential equation is transformed into a system of linear algebraic equations. A closed-form approximation of the fluid pressure is then obtained. The convergence analysis of the proposed method is presented, and its stability is demonstrated numerically using a perturbation approach. The accuracy of the method is validated by comparing the results of the proposed operational matrix scheme with existing methods in the literature. Numerical results indicate that the influence of the fractional-order derivative on pressure distribution depends on both the boundary conditions and the permeability of the medium.</p>