The aim of this thesis was to test the efficiency in practice of an analytical propagator with collision detection for N-body Keplerian systems. This can be used to simulate the evolution of a protoplanetary disk, which gives insight into how planetary systems form. The analytic propagator calculates collisions one by one, while a numerical propagator would compute each time step. The idea of using the analytic propagator is that collisions are rare in astronomical scales, such that jumping from collision to collision and calculating it, is more efficient than calculating all the time steps that are between collisions. Simplifying the orbits of the planetesimals into perfect Keplerian orbits, analytical solutions exist which are used by the analytic propagator. In this thesis, the runtimes of simulations were measured as well as other properties directly related to the runtime. The overall efficiency of the algorithm with respect to N seemed to be O(N3), which is one power less than previously predicted. The prediction was that the runtime of the full simulation is O(N2ε + N4s3/Ia3). Here ε is the maximum eccentricity, s/a is the ratio of a planetesimal's radius to the semi-major axis of its orbit, and I is the maximum inclination. This was calculated by estimating the total number of collisions to be O(N2s2/Ia2) and the runtime for each collision to be O(N2s/a). But the number of collisions turns from quadratic to linear in N, implying that above a certain N almost all planetesimals collide, which reduces the power of N by one. For comparison, the octree code has an algorithmic efficiency of O(N logN) per time step, and the number of steps for a fixed integration time grows as O(N4/3 log N).