The development of an analytically iterative method for solving steady-state as well as unsteady-state problems of cosmic-ray (CR) modulation is proposed. Iterations for obtaining the solutions are constructed for the spherically symmetric form of the CR propagation equation. The main solution of the considered problem consists of the zero-order solution that is obtained during the initial iteration and amendments that may be obtained by subsequent iterations. The finding of the zero-order solution is based on the CR isotropy during propagation in the space, whereas the anisotropy is taken into account when finding the next amendments. To begin with, the method is applied to solve the problem of CR modulation where the diffusion coefficient $\kappa$ and the solar wind speed $u$ are constants with an Local Interstellar Spectra (LIS) spectrum. The solution obtained with two iterations was compared with an analytical solution and with numerical solutions. Finally, solutions that have only one iteration for two problems of CR modulation with $u = constant$ and the same form of LIS spectrum were obtained and tested against numerical solutions. For the first problem, $\kappa$ is proportional to the momentum of the particle $p$, so it has the form $\kappa=k_0\eta$, where $\eta=\frac{p}{m_0c}$. For the second problem, the diffusion coefficient is given in the form $\kappa=k_0\beta\eta$, where $\beta=\frac{v}{c}$ is the particle speed relative to the speed of light. There was a good matching of the obtained solutions with the numerical solutions as well as with the analytical solution for the problem where $\kappa = constant$.