ONE of the most frequently used methods for solving sets of non-linear equat ions g(x) = d, where g(x) is a fairly smooth vector function of a vector argument x with N components g i(x), is Newton's method (see [1, 2]). The computational scheme for this method can be written as: Γ( y ( k) ) z ( k) = g( y ( k) ) − d, y ( k + 1) = y ( k) − z ( k) ( k = 0,1,2,…). Here Γ( y ( k) ) is the Jacobian matrix Γ(y) = ( ∂g ∂y 1 , ∂g ∂y 1 ,2., ∂g ∂y N ) , taken with y = y ( k) . In many cases in practice, in particular in thesolution of non-linear boundary value problems, some other computational scheme is used for the method which replaces the matrix Γ(y (k) in (2) by the difference matrix (R k,y (k=(g(y (k)+h ke 1)−g(y (k)/h k;…;(g(y (k)+h ke N)−g(y (k)/h k where e i is the i-th coordinate unit vector. This scheme takes the form Rh k , y ( k)) z( suk) = g( y ( k) )− d, y ( k+1) = y ( k) − z ( k) (3) ( k = 0,1,2,…), h k are real numbers from some given convergent sequence. In future when we refer to Newton's method we shall mean its modification (3). The convergence of the iterative processes (2) and (3) is generally proved on the assumption that the matrices Γ( y) and R( h, y) are nonsingular in the neighbourhood of the solution x ∗ of (1) (see for instance [1]). In [3] the proof of the convergence of process (2) is given also for the case where the matrix Γ( y) can be singular in the neighbourhood of the point x ∗, but has a constant rank, and all systems in (2) are solvable. In the present paper we shall consider the question of the realization of the iterative process (3) and its convergence in the case where the matrix R( h, y) has a variable rank within the limits 1 to n in the neighbourhood of the solution x ∗, i.e. we shall assume only that the matrix R( h, y) is non-zero in this neighbourhood. In a less general statement this question has been considered for process (2) in [4]. We shall explain what computational difficulties arise in using Newton's method with an unbounded matrix RR −1( h, y) in the neighbourhood of the solution x ∗. Let x ∗ be an isolated singular point of the matrix Γ( y) (i.e. the matrixΓF( x ∗) is singular), and let limR(h k,y (k)) = Γ(x∗) . k→∞ Two cases are then possible. In the first the matrices R( h k , y ( k) ) tend to the singular matrix, remaining non-singular matrices on each iterative process (3). Here we encounter the need to solve with large k systems, with ill-conditioned matrices. The second case consists of the possibility of the appearance of singular matrices R( h k , y ( k) ) on some iterations in (3) (or on all with large k). Consequently in the realization of process (3) in this case we must be able to solve sets of linear equations with a singular matrix, if these systems are solvable, or to reduce the system to solvable form if it has no solution.