  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 $\alpha_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 $5 GeV^2 \leq Q^2 \leq
25 GeV^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 $\lambda_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.
