In this paper, a numerical algorithm is presented to obtain approximate solution of distributed order integro-differential equations. The approximate solution is expressed in the form of a polynomial with unknown coefficients and in place of differential and integral operators, we make use of matrices that we deduce from the shifted Legendre polynomials. To compute the numerical values of the polynomial coefficients, we set up a system of equations that tallies with the number of unknowns, we achieve this goal through the Legendre–Gauss quadrature formula and the collocation technique. The theoretical aspects of the error bound are discussed. Illustrative examples are included to demonstrate the validity and applicability of the method
