We analytically solve the full next-to-leading logarithmic Balitsky-Kovchegov equation in the saturation regime, which includes corrections from quark and gluon loops, and large double transverse logarithms. The analytic result for the $S$-matrix in the saturation regime shows that the linear rapidity decrease with rapidity of the exponent in running coupling case is replaced by rapidity raised to power of 3/2 decreasing. The collinearly-improved Balitsky-Kovchegov equation are also analytically solved in the saturation regime. It shows that the double collinear logarithms do not contribute to the $S$-matrix and the solution is the same as the one obtained from the leading order Balitsky-Kovchegov equation. The numerical solutions to the leading order and full next-to-leading logarithmic Balitsky-Kovchegov equations are performed in order to test the analytic results derived in the saturation regime.