Herein, a spectral Galerkin method for solving the fractional Rayleigh–Stokes problem involving a nonlinear source term is analyzed. Two kinds of basis functions that are related to the shifted sixth-kind Chebyshev polynomials are selected and utilized in the numerical treatment of the problem. Some specific integer and fractional derivative formulas are used to introduce our proposed numerical algorithm. Moreover, the stability and convergence accuracy are derived in detail. As a final validation of our theoretical results, we present a few numerical examples.