Propelled by the discovery of the first exoplanet thirty years ago, the scientific community has rallied and made tremendous strides towards a full understand of the formation and evolution of planetary systems. During this period, over 4, 300 confirmed exoplanets have been detected, and the resulting dataset has driven a revolution by allowing for new formation theories to be proposed and tested. Despite these advances, there is still much about these processes that remains unknown. Numerical n-body simulations of planetary systems are now commonly used to push these frontiers. When performing these investigations, the numerical integration process still poses a very specific set of challenges and therefore demands the continued development of state-of-the-art tools. There are three key novel components to this thesis. Firstly, I perform a detailed analysis of numerical integrators built around multistep collocation methods to quantify their numerical performance over a wide subset of their possible configuration space. Highly favourable performance is observed when specific configurations are applied to globally stiff problems. Secondly, I present my new tool, the Terrestrial Exoplanet Simulator (TES), a novel n-body integration code for the accurate and rapid propagation of planetary systems in the presence of close encounters. TES builds upon the classic Encke method and integrates only the perturbations to Keplerian trajectories to reduce both the error and runtime of simulations. A suite of numerical improvements is presented that together make TES optimal in terms of growth of energy error. Lower runtimes are found in the majority of test problems when compared to direct integration using other leading tools. Finally, using TES, I perform a large simulation campaign to further understand the stability of compact three-planet systems. This work addresses a key limitation in the majority of stability studies by using TES to integrate precisely up to the first collision of planets. Integrations span up to a billion orbits to explore a wide parameter space of initial conditions in both the co-planar and inclined cases. I calculate the probability of collision over time and determine the probability of collision between specific pairs of planets. I find systems that persist for over 108 orbits after an orbital crossing and show how the post-crossing survival time of systems depends upon the initial orbital separation, mutual inclination, planetary radius, and the closest encounter.