This paper revisits the use of the Taylor series method for the numerical integration of ODEs and DAEs. The numerical method is implemented using an efficient variablestep variable-order scheme. Several numerical tests comparing with well-established numerical codes are presented.