In the present study, we introduce two new operational matrices of fractional Legendre function vectors in the sense of generalized Caputo-type fractional derivative and generalized Riemann-Liouville-type fractional integral operators. The derivative and integral operational matrices developed in the sense of Caputo and Riemann-Liouville operators are special cases of our proposed generalized operational matrices for (Formula presented.). Then, we present a numerical method that is dependent on the generalized derivative and integral operational matrices. The applicability and accuracy of the presented method is tested by solving various problems and then comparing the results obtained otherwise by using various numerical methods including spectral collocation methods, spectral Tau method, stochastic approach, and Taylor matrix approach. Moreover, our presented method transforms the problems into Sylvester equations that are easily solvable by using MATLAB or MATHEMATICA. We believe that the newly derived generalized operational matrices and the presented method are expected to be further used to formulate and simulate many generalized Caputo-type fractional models.