Galerkin spectral method. Operational matrices for derivatives have been used and
extended to deal with the system of ordinary differential-algebraic equations. An algorithm
via those matrices has been designed. The accuracy and efficiency of the proposed
algorithm had been shown by two techniques, theoretically, via the boundedness of the
approximated expansion and numerically through numerical examples.