Numerical methods for integrating gravitational 3- and 4-body systems are investigated and tested. The methods employ multiple-pair regularization schemes for N-body systems proposed by Aarseth, Zare, and Heggie, which use the Kustaanheimo-Stiefel transformation for regularizing 2-body collisions, in conjunction with a number of different time transformations between “physical” and “parameter” time. These transformations can be chosen so as to make the singularity in the equations of motion, caused by many-body collisions, as mild as possible. The various time transformations are tested on both 3- and 4-body systems by comparing the numerical with known analytical solutions, and by time reversal of the integrations through many-body close encounters. In addition, the technique of Zare and Szebehely for stabilizing numerical integrations is investigated for each of the time transformations. It is found that when the time transformation is a function of the interparticle separations only, stabilization leads to a significant improvement in performance for those transformations in which the singularity occuring due to many-body collisions is strong (algebraic), but has no decisive effect when the singularity has been softened by an optimal choice of transformation. The most satisfactory time transformation appears to be one involving the Lagrangian, together with the regularization scheme proposed by Aarseth and Zare. Computer programs for binary-single star and binary-binary scattering have been developed and are described. They can be used in an extensive project for determining scattering cross sections with any of the above methods. They are used here to compare the performance of these methods, for a fixed set of intial conditions, on scattering involving “hard” binaries, in which strong resonances can occur. It is found that the outcome of the scattering, for example, the identity of the escaping particle(s), can vary with method, thus reflecting the inherent instability of the N-body problem.