The Taylor formula is used directly in a method of numerical integration of the n-body problem of celestial mechanics; the derivatives in the expansion of the coordinates are calculated successively at each integration step according to the generalized Steffensen rule. The proposed method is the most precise of all numerical methods based on the predetermined part of the Taylor series. The method is used with a variable number of derivatives at each integration step and also with a variable step. The cumulative error in the coordinates increases more slowly in our method than in any other. We can apply the method to the study of the motion of a comet or minor planet, taking into account the perturbations by eight major planets; the method allows for the simultaneous integration of a great number of objects of zero mass.