Abstract
In this paper, an automatic symbolic version of the recurrent power series method has been established for integrating the equations of motion of comets taking into account the perturbations of the giant planets and the non-gravitational effects terms. This new version has the advantage of selecting the coefficients of the power series solutions from the equations of motion of the respective body directly. In this case, it is no need to make derivations of many recurrence relations for the power expansions both of cometary and planetary data. The program is also provided by optimum step size integrations and elaborated using Mathematica software packages. Application to the orbital motion of short-period comets is given.