Third-order aberrations need the paraxial raytracing of a marginal as well as principal ray through the entire system to be determined. When the optical system has GRIN media, even the paraxial raytrace could be cumbersome. This work presents the use of the ABCD matrix for obtaining the raytracing of those rays and, thus, it could ease the calculation of the third-order aberration coefficients, i.e. contribution of GRIN to surface refraction and transfer, for each individual surface of GRIN optical elements. A procedure to be implemented in any modern commercial optical software has been presented and evaluated for simple systems.