We present an implementation of the Wisdom-Holman integrator for simulating gravitational dynamics in planetary systems: systems with one dominant central mass and \(N\) orbiting bodies, such as the Solar System. The Wisdom-Holman integrator models the motion of non-central bodies as unperturbed Kepler orbits and integrates gravitational interactions between orbiting planets as weak perturbations. Two methods to advance a body along its orbits are investigated: one using coordinate transformations (Method A) and one based on the \(f\) and \(g\) functions (Method B). Method B is shown to be significantly faster than Method A, without a significant loss of accuracy, making it the preferred method for most simulations. Simulations of the Solar System using large time steps, including time steps exceeding the orbital period of some planets, are explored to determine whether the increase in computational speed justifies the loss in accuracy. Simulations with a fixed Sun and with a dynamic Sun are considered. Results indicate that for fixed Sun simulations, while accurate simulations require small time steps, larger steps still capture the qualitative behaviour of the system. The step size of simulations with a dynamic Sun is limited by \(\Dt\st{max} = T_\mercury / 6\). However, for small step sizes, dynamic Sun simulations accurately describe the Solar System and the restricted three-body problem in which resonance occurs. This is achieved without the use of Jacobian coordinates, which are commonly used in implementations of the Wisdom-Holman integrator. The fixed Sun and dynamic Sun simulations are shown to conserve energy, suggesting an accurate description of the system simulated.