In this work we have solved the nonlinear GLR-MQ evolution equation upto next-to-leading order (NLO) by considering NLO terms of the gluon-gluon splitting functions and running coupling constant α s (Q 2 ). Here, we have incorporated a Regge-like behaviour of gluon distribution in order to obtain a solution of the GLR-MQ equation in the range of 5GeV 2 ≤ Q 2 ≤ 25GeV 2 . We have studied the Q 2 evolution of the gluon distribution function G(x, Q 2 ) and its nonlinear effects at small-x. It can be observed from our analysis that the nonlinearities increase with decrease in the correlation radius R of two interacting gluons, as expected. We have compared our result of G(x, Q 2 ) as Q 2 increases and x decreases, for two different values of R, viz. R= 2 GeV −1 and 5 GeV −1 . We have also checked the sensitivity of the Regge intercept λ G on our results. We compare our computed results with those obtained by the global analysis to parton distribution functions (PDFs) by various collaborations where LHC data have been included viz. PDF4LHC15, NNPDF3.0, ABM12 and CT14. Besides we have also shown comparison of our results with HERA PDF data viz. HERAPDF15.