In this work 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 and are part of the Tau Toolbox - a Matlab library for the Lanczos' tau method. 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.