subroutine Newton(SUB, x0, eps, lnit, rad, kod, n_div)
! Solutia ecuatiei F(x)=0, prin metoda Newton.
! Parametri:
!	SUB	:numele subrutinei care calculeaza functia f =f(x) si derivata "f1 =f'(x)".
!		 Se declara EXTERNAL in programul care apeleaza NEWTON.
!	X0	:aproximatia initiala
!	eps	:toleranta
!	lnit: -Intrare =nr.maxim de iteratii; -Iesire =nr.efectiv de iteratii.
!	rad = solutia ecuatiei
!	kod	: -INTRARE: cod de tiparire a iteratiilor: 0 = Nu; #0 = Da.
!          kod >=2: se tipareste si corectia curenta.
!       : -IESIRE : cod de terminare/eroare:
!		 0:	toleranta "eps" atinsa;
!		 1:	depasirea nr.maxim de iteratii admis "lnit";
!        2: Proces stationar;
!		-1: divergenta
!       -2:	derivata functiei in X0 este ~ nula.
!
!* Exemplu de subrutina SUB:
! subroutine SUB(x,f,f1)
! common p	! optional, in cazul cand F include un parametru "p".
!			! "p" se citeste in programul apelant.
! f =...
! f1=...
! end

logical print_it, print_cor
logical stationar
integer width

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

logical print_div


print_it = kod /=0
print_cor = kod >=2
print_div =.false.

if(print_it) then
print *
    if(.not.print_cor) then 
    write(1, '(a)') ' Iteratia     x'
    else
    write(1, '(a)') ' Iteratia     x               Corectia curenta'
    endif
endif

! Initialize:
aa =0.
aa(0) =x0

iter =0

x =x0
dif_ant =100.
!n_div =0

100	iter =iter +1
call sub(x,f,f1)

! Test derivata nula:
if(abs(f1) < 1.E-10) then		! tiny(f1)
print *, 'f-derivat =', f1
write(1,*) 'f-derivat =', f1

kod=-2
goto 500
endif

x  =x0 -f/f1		!dif

dif =x-x0

iwi =width(iter)
if(print_it) then
	if(.not.print_cor) then 
    print 90, iter, x
	write(1,90) iter, x
	
	else
    print 91, iter, x, dif
	write(1,91) iter, x, dif
	endif
endif

90	format(/, I9, t13,1pG14.7,\)
91	format(/, I9, t13,1pG14.7, t30,1pG14.7,\)

! EPS = atins:
IF(abs(dif) <= eps) THEN
kod =0
goto 500
ENDIF

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

IF(per >0) 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

! Divergenta locala:
IF(abs(dif) > abs(dif_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 numar de iteratii:
IF(iter < lnit) THEN
goto 200

ELSE
	if(abs(dif) > abs(dif_ant)) then
	kod =-1		! Divergenta
	else
	kod=1		! LNIT depasit
	endif

goto 500
ENDIF

! Noua iteratie
200 continue
dif_ant =dif
x0 =x
goto 100


500 continue
rad =x
lnit =iter

if(.not. print_it .and. kod ==2 .and. per <=50) then
write(1,*)
do ii =i, iter
write(1, '(I6, t14, G16.7)') ii, aa(ii)
enddo
endif

!print '(/, a, 1pg14.7)','  Diferenta curenta = ', dif
if(.not.print_cor) then 
print '(a, 1pg14.7)','  Corectia curenta = ', dif
endif

!if(print_it .or. kod ==2) write(1,'(/)')
if(.not.print_cor) write(1, '(a, 1pg14.7)') '  Corectia curenta =', dif


END
