The long-term simulation of planetary systems poses significant challenges due to the inherently chaotic and non-integrable nature of gravitational interactions in the N-body problem. This thesis examines the Wisdom–Holman symplectic integration scheme, a method specifically designed for nearly integrable systems. This scheme separates the dominant Keplerian motion from weaker perturbative forces, enabling stable integration over astronomical timescales. Emphasis is placed on understanding the practical limitations and capabilities of this method when using large time steps, particularly in the presence of mean-motion resonances and step-size resonances. Through extensive numerical experiments using the Rebound simulation package, the scaling behavior of integration errors is characterized, revealing a transition from secondorder to lower-order error regimes at large time steps. This shows that in certain systems, time steps significantly larger than the shortest orbital period can still yield acceptable accuracy. However, in systems with mean-motion resonance, strong sensitivity to step-size resonances is observed, requiring careful step-size selection. A comparison between Jacobi and Democratic Heliocentric coordinates shows that the former performs best when orbits are nested, while the latter is better suited to systems with crossing or unordered orbits. These findings provide practical guidelines for applying Wisdom–Holman integration effectively across a range of dynamical regimes.