The paper is concerned with iteratively solving systems of nonlinear equations $g(x)=0$, where $g$ is a continuously differentiable mapping in $n$-dimensional real space. It is supposed that the systems are large-scale systems for which the Jacobian is not available or requires a prohibitive amount of storage. The author extends the conjugate gradient method to solve the system (a problem equivalent to an unconstrained optimization-minimization problem). For this, the known Polak-Ribiere-Polyak conjugate gradient direction , as a new line search direction, is used [see {\it E. Polak} and {\it G. Ribière}, Rev. Franç. Inform. Rech. Opér. 3, No. 16, 35--43 (1969;

Zbl 0174.48001) and {\it B. T. Polyak}, U.S.S.R. Comput. Math. Math. Phys. 9(1969), No. 4, 94--112 (1971); translation from Zh. Vychisl. Mat. Mat. Fiz. 9, 807--821 (1969;

Zbl 0229.49023)]. The author proposes the algorihm DFCGNE (Derivative Free Conjugate Gradient for Nonlinear Equations) for solving nonlinear systems and also, modification of this algorithm, called M-DFCGNE method, in the case of nonmonotone objective functions. Under some reasonable conditions, the global convergence of these algorithms is proved. Numerical experiments and comparisons with other methods are discussed.