Emulating the long-term evolution of N-body gravitational systems is a significant challenge for standard machine learning models, which typically fail to respect fundamental conservation laws, leading to unphysical and unstable trajectory predictions. We address this by developing a Symplectic Neural Ordinary Differential Equation framework designed to learn the underlying conservative vector field governing the dynamics. Our model parameterizes the system's Hamiltonian using a permutation-invariant graph neural network, from which forces are derived via automatic differentiation to ensure they are curl-free. Crucially, we embed a differentiable leapfrog integrator directly into the training loop, which constrains the learned dynamics to be symplectic. Training is performed on trajectory snapshots from simulations of 50-particle virialized Plummer spheres, where a gravitational softening length is incorporated as a fixed physical prior and a curriculum learning strategy is employed to handle the system's multi-scale density. This approach transforms the learning problem from brittle state-to-state regression into the robust emulation of a continuous Hamiltonian flow. By construction, the learned dynamics preserve the geometric structure of the phase space, exhibiting long-term energy stability, time-reversibility, and phase-space volume conservation. The resulting emulator generalizes to systems with different particle counts, demonstrating that explicitly encoding physical symmetries is a more effective path to building robust models for chaotic physical systems than purely minimizing trajectory error.