This research paper focused on the solution of systems of fractional integro-differential equations (FIDEs) of the Volterra type with variable coefficients. The proposed approach combined the tau method and shifted Gegenbauer polynomials in a matrix form. The investigation of the existence and uniqueness of solutions for these systems was carried out using Krasnoselskii's fixed point theorem. The equations employed Caputo-style derivative operators, and to minimize computational operations involving derivatives and multiplications, integral and product operational matrices were derived. By introducing suitable polynomial approximations and employing the tau spectral method, the original system of FIDE was transformed into an algebraic system. Solving this algebraic system provided approximate solutions to the main system. Error bounds were computed in the Gegenbauer-weighted Sobolev space. The proposed algorithm was implemented and tested on two systems of integro-fractional differential equations to demonstrate its efficiency and simplicity. By varying the parameter $\sigma$ in the Gegenbauer polynomials, the impact of this variation on the approximate solutions can be observed. A comparison with another method utilizing the block-by-block approach was also presented.