In this paper, numerical technique based on monic Laguerre polynomials is proposed to obtain approximate solutions of initial value problems for systems of fractional order integro-differential equations (FIDEs). Operational fractional integral matrix is constructed. This operational matrix is applied together with the monic Laguerre Tau method to solve systems of FIDEs. This systems of FIDEs will be transformed into a system of algebraic equations which can be solved easily. Numerical results and comparisons with other methods are also presented to show the efficiency and applicability of the proposed method.