<p>In this manuscript, we develop an efficient numerical scheme for solving time distributed-order fractional integro-differential equations. The proposed method utilizes orthogonal generating polynomials (OGPs) and Chebyshev-polynomial-based operational matrices. First, we construct OGPs associated with a distributed-order weight function defined in the Caputo sense. Using these polynomials, we derive operational matrices corresponding to both the integral and distributed-order fractional derivative terms, incorporating the Legendre–Gauss quadrature rule for numerical integration. The original complex problem is then transformed into a system of linear algebraic equations using the standard tau method, combined with uniformly distributed random number as collocation points to enforce the initial and boundary conditions. The resulting approximate solutions demonstrate the accuracy and efficiency of the proposed approach. Furthermore, we establish error bounds, conduct a detailed convergence analysis, and verify the numerical stability of the scheme. Several illustrative examples are provided to highlight the reliability and computational performance of the method. For comparison, results obtained using the Chebyshev polynomial based operational matrix method are also presented.</p>