The aim of this paper is to encourage the use of orbital integrators in the classroom to discover
and understand the long term dynamical evolution of systems of orbiting bodies. We show how to
perform numerical simulations and how to handle output data in order to reveal the dynamical
mechanisms that dominate the evolution of arbitrary planetary systems in timescales of millions of
years using a simple but efficient numerical integrator. Through some examples we reveal the
fundamental properties of planetary systems: the time evolution of the orbital elements, the free
and forced modes that drive oscillations in eccentricity and inclination, the fundamental
frequencies of the system, the role of the angular momenta, the invariable plane, orbital
resonances, and the Kozai–Lidov mechanism.