In this paper we proposed the homotopy approach for solving the nonlinear Balitsky-Kovchegov evolution equation with running QCD coupling. The approach consists of two steps. First, is the analytic solution to the nonlinear evolution equation for the simplified, leading twist kernel. Second, is the iteration procedure that allow us to calculate corrections analytically or seminumerically. For the leading twist kernel it is shown that the first iteration leads to <= 1% accuracy. The zeta 1/4 - 4Nc dipole size, Qs is the saturation scale) and geometric scaling behavior of the scattering amplitude are discussed as well as the dependence on the value of the infrared cutoff.