We calculate the NLO corrections for the gluon fragmentation functions to a heavy quark-antiquark pair in ${^{1}hspace{-0.6mm}S_{0}^{[1]}}$ or ${^{1}hspace{-0.6mm}S_{0}^{[8]}}$ state within NRQCD factorization. We use integration-by-parts reduction to reduce the original expression to simpler master integrals (MIs), and then set up differential equations for these MIs. After calculating the boundary conditions, MIs can be obtained by solving the differential equations numerically. Our results are expressed in terms of asymptotic expansions at singular points of $z$ (light-cone momentum fraction carried by the quark-antiquark pair), which can not only give FFs results with very high precision at any value of $z$, but also provide fully analytical structure at these singularities. We find that the NLO corrections are significant, with K-factors larger than 2 in most regions. The NLO corrections may have important impact on heavy quarkonia (e.g. $eta_c$ and $J/psi$) production at the LHC.