The time fractional subdiffusion equation (FSDE) as a class of anomalous diffusive systems has obtained by replacing the time derivative in ordinary diffusion by a fractional derivative of order $0<\alpha <1$. Since analytically solving this problem is often impossible, proposing numerical methods for its solution has practical importance. In this paper, an efficient and accurate Galerkin method based on the Legendre wavelets (LWs) is proposed for solving this equation. The time fractional derivatives are described in the Riemann–Liouville sense. To do this, we first transform the original subdiffusion problem into an equivalent problem with fractional derivatives in the Caputo sense. The LWs and their fractional operational matrix (FOM) of integration together with the Galerkin method are used to transform the problem under consideration into the corresponding linear system of algebraic equations, which can be simply solved to achieve the solution of the problem. The proposed method is very convenient for solving such problems, since the initial and boundary conditions are taken into account, automatically. Furthermore, the efficiency of the proposed method is shown for some concrete examples. The results reveal that the proposed method is very accurate and efficient.