  This paper provides an efficient recursive approach of the spectral Tau
method to approximate the solution of system of generalized Abel-Volterra
integral equations. In this regards, we first investigate the existence,
uniqueness as well as smoothness of the solutions under assumption on the given
data. Next, from a numerical perspective, we express approximated solution as a
linear combination of suitable canonical polynomials which are constructed by
an easy to use recursive formula. Mostly, the unknown parameters are calculated
by solving a low dimensional algebraic systems independent of degree of
approximation which prevent from high computational costs. Obviously, due to
singular behavior of the exact solutions, using classical polynomials to
construct canonical polynomials, leads to low accuracy results. In this
regards, we develop a new fractional order canonical polynomials using
M\"untz-Legendre polynomials which have a same asymptotic behavior with the
solution of underlying problem. The convergence analysis is discussed, and the
familiar spectral accuracy is achieved in $L^{\infty}$ norm. Finally, the
reliability of the method is evaluated using various problems.
