—in this paper, present a computational method for solving Fredholm integral equations of the second kind. The method based on the application of the shifted Legendre polynomials in matrix forms. We create a technique for extracting the Legendre coefficients of each polynomial away so that each Legendre polynomial is rewritten in the form of its coefficient's matrix multiplied by the monomial basis function matrix. This technique significantly reduces the round off errors. By using this technique, the unknown and the data functions are expressed in two forms; each consists of three matrices. The kernel is approximated twice relevant to its two variables so that it is transformed into a form consists of a product of five matrices. By substituting by the approximate unknown function into the left and the right sides of the integral equation, we obtain an algebraic linear system of the equations without applying the collocation points. Moreover, we adapted the Gauss–Quadrature rule in an adjustment form and applied it for computing the resulted integrals. The convergence in the mean of the approximate solution and the kernel are proved. Additionally, the maximum norm error is studied, and it is found equal to zero. Numerical results are obtained for five examples to clarify the simplicity, efficiency, and reliability of the method. The obtained solutions are equal or rapidly converge to the exact solutions.