Ray tracing in gradient-index (GRIN) media has been traditionally performed either by using the analytical or numerical solutions to the Eikonal equation or by creating a layered medium where Snell's law is calculated in each layer. In this paper, an exact general method to perform ray tracing in GRIN media is presented based on the invariants of the system as stated by Fermat's principle when the media presents symmetries. Its advantage, compared with other methods reported in the literature, relies on its easy implementation. Besides the GRIN distribution and the initial conditions of the incident ray, once the invariants of the system are stated the resulting math is simple to solve and interpret. To benchmark the algorithm, ray tracing in typical cases of GRIN media is calculated, finding minimal discrepancies between the analytical solutions and our simulations. The used media are axial refractive index and parabolic index fiber and lenses with spherical gradient-index symmetry, such as: Luneburg's, Gutman's, generalized Maxwell's Fish-eye, Eaton's, and concentrator lenses. Our method can be further applied to distributions with symmetries associated with other common curvilinear orthogonal coordinate systems, in particular to those associated to the separability of the Helmholtz equation that would allow us to investigate wave optics in these GRIN media with the associated geometries.