In this chapter, a novel numerical method based on triangular orthogonal functions is proposed for the numerical solutions of fractional-order integro-differential equations such as Fredholm integro-differential equations of fractional order, Volterra integro-differential equations of fractional order, and Fredholm-Volterra integro-differential equations of fractional order. It is theoretically shown that there exists a unique solution to the general form of the system of fractional-order integro-differential equations considered in this chapter. Convergence analysis is conducted to prove that in the limit of step size tends to zero, the proposed numerical method can converge the approximate solution to the exact solution of system of fractional integro-differential equations considered. Numerical examples as well as physical process models involving fractional-order integro-differential equations are solved to demonstrate the effectiveness of the proposed numerical method.