In this work, we discuss an operational matrix approach for introducing an approximate solution of the fractional subdiffusion equation (FSDE) with both Dirichlet boundary conditions (DBCs) and Neumann boundary conditions (NBCs). We propose a spectral method in both temporal and spatial discretizations for this equation. Our approach is based on the space-time shifted Legendre tau-spectral method combined with the operational matrix of fractional integrals, described in the Riemann–Liouville sense. The main characteristic behind this approach is to reduce such problems to those of solving systems of algebraic equations in the unknown expansion coefficients of the sought-for spectral approximations. In addition, this approach is also investigated for solving the FSDE with the variable coefficients and the fractional reaction subdiffusion equation (FRSDE). For conforming the validity and accuracy of the numerical scheme proposed, four numerical examples with their approximate solutions are presented. Also, comparisons between our numerical results and those obtained by compact finite difference method (CFDM), Box-type scheme (B-TS), and FDM with Fourier analysis (FA) are introduced.