A mathematically coherent approach to the stability of orbits in astrophysical problems, that preserves the energy and the symplectic structure of the equations of motion is called Hamiltonian flow analysis. A Hamiltonian-based model of multi-body gravitational interactions, based on the canonical equations of motion, and symplectic integration, is presented to study the long-time evolution of the orbits under perturbative forces. The technique of “phase-space trajectory”, Lyapunov stability analysis and Poincaré section are used to distinguish between stable and chaotic behaviours in orbit. With numerical experiments, the accuracy of the orbit prediction is shown to be 96.8%, Hamiltonian energy conservation is 98.14% and stability in phase space is 94.7% with a 17.9% reduction in the cumulative trajectory error as compared to conventional Runge–Kutta integration. The data shows that there are improvements to the numerical stability and physically consistent long-term orbital predictions. A mathematical model for celestial mechanics, a reliable model for planetary system analysis and spacecraft trajectory design is proposed.