  In a Keplerian system, a large number of bodies orbit a central mass.
Accretion disks, protoplanetary disks, asteroid belts, and planetary rings are
examples. Simulations of these systems require algorithms that are
computationally efficient. The inclusion of collisions in the simulations is
challenging but important. We intend to calculate the time of collision of two
astronomical bodies in intersecting Kepler orbits as a function of the orbital
elements. The aim is to use the solution in an analytic propagator ($N$-body
simulation) that jumps from one collision event to the next. We outline an
algorithm that maintains a list of possible collision pairs ordered
chronologically. At each step (the soonest event on the list), only the
particles created in the collision can cause new collision possibilities. We
estimate the collision rate, the length of the list, and the average change in
this length at an event, and study the efficiency of the method used. We find
that the collision-time problem is equivalent to finding the grid point between
two parallel lines that is closest to the origin. The solution is based on the
continued fraction of the ratio of orbital periods. Due to the large jumps in
time, the algorithm can beat tree codes (octree and $k$-d tree codes can
efficiently detect collisions) for specific systems such as the Solar System
with $N<10^8$. However, the gravitational interactions between particles can
only be treated as gravitational scattering or as a secular perturbation, at
the cost of reducing the time-step or at the cost of accuracy. While
simulations of this size with high-fidelity propagators can already span vast
timescales, the high efficiency of the collision detection allows many runs
from one initial state or a large sample set, so that one can study statistics.
