Summary: The universal formulation for the perturbed two-body problem is generalized to cover all gravitational \(N\)-body problems involving a dominant central mass. Its efficiency, when compared to conventional numerical integration, is shown in several examples. The convergence and numerical stability of the method is discussed, and a universal state transition matrix is obtained, which can be used either in a process of differential correction of an orbit or, as in the present case, to obtain an accurate estimation of global errors.