subroutine Newton_Sys (n, x, lnit, xtol, f, kod, n_div)
! Solutia unui sistem de "n" ecuatii neliniare: F(X) =0.
! Metoda NEWTON - cu Jacobianul analitic
! fcn, jfcn: subrutne care calculeaza F(X) si jacobianul JF(X)
! 
! Parametri:
! n: numarul ecuatiilor
! lnit: INTRARE/IESIRE: numarul limita/ numarul efectiv de iteratii
! x: tabloul necunoscutelor. La intrare contine aproximatiile initiale.
! f: vectorul coordonatelor functirei F(X)
! xtol: toleranta la delta_x
! kod: INTRARE: cod de tiparire iteratii: =0/#0: Da/Nu.
!   IESIRE: cod de incheiere a iteratiei:
!		0: XTOL atins
!		1: LNIT depasit
!		2: Proces stationar;
!		-1: divergenta
!       -2:	Pivot < prag in eliminarea Gauss (Jacobian singular)
!

real(4) X(N), F(N), JF(N,N), w(n)	! w = corrections to X.
real(4) A(n,n), b(n)
LOGICAL print_it
real norm$

real prag/1E-7/		! pt. Elim

real xall(n, 0:lnit)	! Tablou pentru toate val. - pt. proces stationar
integer period, per

logical print_div

integer width
character(22) restart_text(2)

!epsX(x1) =1./2**(24 -exponent(x1))

print_it =(kod /= 0)
print_div =.false.

restart_text(1) ='Repornire iteratie:'
restart_text(2) ='Restarting Iteration:'
ilang =1

! Iteratie:
iter =0
xall =0.
xall(:,0) =x
wnorm_ant =100.
n_div =0

iRestart =0
100 iter =iter +1
! Calculeaza F si JF:
call fcn(n, x, f)
call jfcn(n, x, jf)

! Matricea A pentru Gauss: JF+ termenii liberi:
A =jf
b =-F(:)

nb =1
kodPrint =0

call ELIM (n, nb, A, b, prag, kodElim, kodPrint)

! Test Jacobian singular:
if(kodElim ==-1) then
kod =-2

!lnit =iter
print*, '** Matrice Jacobian singulara'
!return

goto 500
endif

! Corecteaza X:
w = b
x = x + w

! Tiparire valori curente
	if(print_it) then
	call fcn(n, x, f)		! Pentru F(X)
	if(iter ==1) print 90
	print 91, iter, x, f

	if(iter ==1) write(1,90)
	write(1,91) iter, x, f
	endif

90 format(' Iteratia:', t14, 'x ', t<n*16+13>, 'f')
91 format(i6, t10, <n>G16.7, t<n*16+10>, <n>G16.7)

! Test XTOL:
wnorm =norm$(n, w)

IF(wnorm <= XTOL) THEN
kod =0
!lnit =iter
!if(.not. print_it) call fcn(n, x ,f)
! Tiparire valori curente
	print*,' Norma corectiei:', wnorm
	write(1,*) ' Norma corectiei:', wnorm

if(.not. print_it) then
call fcn(n, x ,f)			! Pentru F(X)
print *, iter, x, f
endif

goto 500
ENDIF

! Test Stationar:
xall(:, iter) =x		! pt. stationar
per =period(n, iter, xall, i)

!IF(per >0) then
IF(per >1) then
iwiter =width(iter)
iwi =width(i)
iwper =width(per)

!	write(*, '(/, a, 1pG14.7, \)') ' Radacina', x
	print '(/, a, I<iwiter>, a, I<iwi>, a)', ' x(', iter, ') = x(', i, ')'
	print '(a, I<iwper>)', ' ** Proces stationar de perioada ', per

!	write(1, '(//, a, 1pG14.7, \)') ' Radacina', x
	write(1, '(/, a, I<iwiter>, a1, \)') ' Iteratia ', iter,': '
	write(1, '(a, I<iwiter>, a, I<iwi>, a)') ' x(', iter, ') = x(', i, ')'
	write(1, '(a, I<iwper>)') ' ** Proces stationar de perioada = ', per

kod =2
goto 500
ENDIF

! Test Divergenta locala:
IF(abs(wnorm) > abs(wnorm_ant)) THEN
!kod =-1 
n_div =n_div+1

if(print_it) then
if(.not. print_div) print *, ' ** Divergenta locala la iteratia:'
print_div =.true.

print *, iter

write(1, '(a)') ' ** Divergenta locala'
endif
ENDIF

! Test nr. de iteratii:
IF(iter == lnit) THEN
kod =1
	print*,' Norma corectiei:', wnorm
	write(1,*) ' Norma corectiei:', wnorm

if(.not. print_it) then
call fcn(n, x, f)
print *, iter, x, f
endif

ELSE	! Noua iteratie

wnorm_ant =wnorm
goto 100

ENDIF


500 lnit =iter

END subroutine Newton_Sys

