This study proposed a method for calculating any arbitrarily linear polarization paraxial optical field propagation in Gradient-index (GRIN) lens. The proposed method uses multiple thin-phase sheets to approximate a GRIN lens. This study also compares the proposed method with the other well-known method. That is, using a single lens equivalent optical system of Fractional Fourier transform (FrFT) to simulate a GRIN lens. The evolution in GRIN lenses of many special beams have been calculated by the FRFT-method. This study use both methods to calculate the Helmholtz-Gauss beam evolution in GRIN lenses of a small and a high gradient constant, respectively. Numerical results shows that the differences between two calculation methods appeared while the GRIN lens of a high gradient constant. This study provides an alternative approach could calculate the linearly polarized field evolution in GRIN lenses with higher precision, which will be useful to the optical design of GRIN lens systems.