  In this paper, we explore the chaotic signatures of the geodesic dynamics for
particles moving in the slowly rotating Hartle-Thorne spacetime; an approximate
solution of vacuum Einstein field equations describing the exterior of a
massive, deformed, and slowly rotating compact object. We employ the numerical
study to examine the geodesics of prolate and oblate deformations for generic
orbits and find the plateaus of the rotation curve, which are associated with
the existence of Birkhoff islands in the Poincare surface of the section, where
the ratio of the radial and polar frequency of geodesics remains constant
throughout the island. We investigate various phase space structures, including
hyperbolic points and chaotic regions in the neighborhood of resonant islands.
Moreover, chaotic behavior is observed to be governed by the stickiness
phenomenon, where chaotic orbits remain attached to stable ones for an extended
duration before eventually diverging and are attracted toward the surface of
the neutron star. The precision of the numerical integration used to simulate
the particle's trajectories plays a crucial role in the structures of the
Poincare surface of the section. We present a comparison of several efficient
structure-preserving numerical schemes of order four applied to the considered
non-integrable dynamical system and we investigate which schemes possess the
canonical property of the Hamiltonian flow. Among the class of non-symplectic
integrators, we employ the explicit Runge-Kutta method and explicit general
linear method with a standard projection technique to project the numerical
solution onto the desired manifold. The projection scheme admits the
integration without any drift from the desired manifold and is computationally
cost-effective. We are concerned with two crucial aspects -- long-term
behaviour and CPU time consumption.
