  In this paper an extension of the spectral Lanczos' tau method to systems of
nonlinear integro-differential equations is proposed. This extension includes
(i) linearization coefficients of orthogonal polynomials products issued from
nonlinear terms and (ii) recursive relations to implement matrix inversion
whenever a polynomial change of basis is required and (iii) orthogonal
polynomial evaluations directly on the orthogonal basis. All these improvements
ensure numerical stability and accuracy in the approximate solution. Exposed in
detail, this novel approach is able to significantly outperform numerical
approximations with other methods as well as different tau implementations.
Numerical results on a set of problems illustrate the impact of the
mathematical techniques introduced.
