subroutine Newton(SUB, x0, eps, lnit, rad, kod)
implicit real(8) (a-h, o-z)

! Solutia ecuatiei F(x)=0, prin metoda Newton.
! Parametri:
!	SUB	:numele subrutinei care calculeaza valoarea lui F si a derivatei F1.
!		 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.
!       : -IESIRE : cod de terminare/eroare:
!		 0:	toleranta "eps" atinsa;
!		 1:	depasirea nr.maxim de iteratii admis "lnit";
!        2: Proces stationar: x2 =x0;
!		-1: divergenta
!       -2:	derivata functiei in X0 este ~ nula.
!
!* Exemplu de subrutina SUB:
! subroutine SUB(f,f1,x)
! common p	! optional, in cazul cand F include un parametru "p".
!			! "p" se citeste in programul apelant.
! f =...
! f1=...
! end

logical print_it
logical stationar

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

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

print_it = kod /=0
if(print_it) print *

! Initialize:
iter =0

x =x0
dif_ant =100.

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

! Test derivata nula:
if(abs(f1) < 1.e-15) then
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
	print 90, iter,x
	write(1,90) iter, x
	endif

90	format(/, I9, t13, 1pG24.15, \)

! EPS = atins:
!IF(abs(dif) <= eps) THEN
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)

	print '(/, a, I<iwiter>, a, I<iwi>, a)', ' x(', iter, ') = x(', i, ')'
	print '(a, I<iwper>)', ' ** Proces stationar de perioada ', per

	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 =-2 

if(print_it) then
print *, ' ** Divergenta locala'
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, G25.16)') ii, aa(ii)
enddo
endif

print '(/, a, 1pg25.16)','  Diferenta curenta = ', dif

if(print_it .or. kod ==2) write(1,'(/)')
write(1, '(a, 1pg25.16)') '  Diferenta curenta =', dif

END
