! y =x -f(x)/f'(x)
! z =y -f(x)/(f(x) -2*f(y))*f(y)/f'(x)
! w =z -f(z)*f[x,y]/div_dif(f,x,z)/div_dif(f,y,z)*(1+ f(z)/f(x))
! x1 =w -f(w)/( div_dif(f,x,w) +( div_dif_2(f,y,x,z) -div_dif_2(f,y,x,w) -div_dif_2(f,z,x,w) )*(x-w) )
!
! Ordin 15:

subroutine Ordin_15(x0, eps, lnit, rad, kod, n_div)
use work
! Solutia ecuatiei F(x)=0, prin metoda de Ordin 15.
! Parametri:
!	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 tiparesc si (marimile/pas): y,z,w si f(y),f(z),f(w)si numitorii la z, w.
!   kod (Iesire): cod de terminare/eroare:
!		 0:	toleranta "eps" atinsa;
!		 1:	depasirea nr.maxim de iteratii admis "lnit";
!        2: |f(x)| <=eps_f
!        3: Proces stationar;
!		-1: divergenta
!       -2:	Numitor < eps_num
!       -3,-4,-5: y = x; z = y; w = z;
!
logical print_it, print_pas
logical stationar
integer width

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

logical print_div
real num    ! numitorul in formule

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

print_it = kod /=0
print_pas = kod >=2

print_div =.false.

if(print_it) then
print *
write(1, '(a)') ' Iteratia    x'
endif


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

iter =0

x =x0
dif_ant =100.

!n_div =0   ! in Main
!eps_num =1E-12

100	iter =iter +1
x0 =x
call Pas_Ordin_15(x0, x, iter,  print_pas, kod)

IF(kod ==-10) THEN
rad =x
RETURN
ENDIF
!x =w -fw/( div_dif(x,w) +( div_dif_2(y,x,z) -div_dif_2(y,x,w) -div_dif_2(z,x,w) )*(x-w) )
!x =w -dif_w

dif = x -x0

iwi =width(iter)
	if(print_it) then
	print 90, iter,x
	write(1,90) iter, x
	endif

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

! EPS = atins:
IF(abs(dif) <= eps) THEN
!
write(1,*) 'dif:', dif
kod =0
goto 500
ENDIF

IF(abs(f(x)) <= eps_f) THEN
kod =2
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 =3
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

if(kod /=2) then
print '(/, a, 1pg14.7)','  Diferenta curenta =', dif

if(print_it) write(1,'(/)')
write(1, '(a, 1pg14.7)') '  Diferenta curenta =', dif
endif

END
