!The name of this pragrom is 3cn.f90
!f(x)=2*x*cos(2*x)-(x-2)**2
!USE External Module interval_arithmetic
!USE External Procedure Simini
!Use Newton Interval_Arithmetic to Solve Page 9 -Problem 3c
! 
program problem_3c_with_interval
   use interval_arithmetic
   implicit none
! Interval [2.0,4.0], 
! Epsilon is error bound 

   real, parameter :: Tv1=2.0, Tv2=4.0, Epsilon=0.00001
   type(interval) :: x, value,x0

   real:: tempf,tempx, in1,in2,Midp,Lmidp,W,Step=0.1
   integer :: i,j ,nn
   real:: NumInv ,Solu
   dimension :: NumInv(10,2),Solu(10)
   
   call simini
   in1=Tv1; in2=Tv1+Step

! This loop is to shrink the interval of solution of Df(x)=0
   i=1; nn=0
   do  
    x=interval(in1,in2)
!	value= 2*cos(2*x)-4*x*sin(2*x)-2*(x-2)   
	value=Df1(x)
	write(*,*) "I=:", i,value
	 if(value%upper.GE.0 .AND. value%lower.LE.0)then   
	  nn=nn+1 
	  NumInv(nn,1)=in1; NumInv(nn,2)=in2
	 endif
	  
     i=i+1
     in1=in2
	 in2=in2+Step
	 
	if (in2.GT.(Tv2+Step).and. nn.EQ.0)then
	 
	  write(*,*) "No sultion of f(x)=0"
	  exit
	endif
	if (in2.GT.(Tv2+Step).and. nn.GT.0)	then
	  Write(*,*)"There is",nn," solutions"
	exit         
	endif
	end do	 

Do j=1,nn
 Print *,"NN=", nn, NumInv(j,1),NumInv(j,2)	 
END DO

!This loop is to find solution of Df(x)=0
 
 Do J=1,nn !1
   Lmidp=(NumInv(1,1)+NumInv(1,2))/2 
   i=0
  DO  !2				
   Write(*,*) "Enter Do2"    
    i=i+1
    x=interval(NumInv(nn,1),NumInv(nn,2))
    Midp=(NumInv(nn,2)+NumInv(nn,1))/2
    W=abs(NumInv(nn,2)-NumInv(nn,1))
    if(abs(Df(Midp)-Df(Lmidp)).LT.Epsilon .OR. W.LT.Epsilon)then 
      Solu(nn)=Midp
      write(*,*) "The solution of Df(x)=0 is ",Midp
	  exit
     else
      x0=interval(Midp-Epsilon,Midp+Epsilon)
	   
	  value=X0-Df1(X0)/DDf(x)
	  
	  x=value.IS.x
	  print *, "Value=", value,  x
 
	  NumInv(nn,1)=x%upper
      NumInv(nn,1)=x%lower
	  Lmidp=Midp
    endif
	if (i.GT.10) then 
	exit
	endif
   	 
   end do ! 2
   write (*,*) "Out D02"
  END DO !   1


   write(*,*) "RESULT"
   
   DO I=1, nn
    Write (*,*) "X=", Solu(i), "     F(x)=", f(Solu(i))
   End DO 

  
   
  If (abs(f(Tv1)).GT.Abs(f(Tv2))) then
   Tempf=abs(f(Tv1))
   Tempx=Tv1
  else
   Tempf=abs(f(Tv2))
   Tempx=Tv1
  endif

 Do i=1,nn
  If (tempf .LT.abs(f(Solu(i)))) then
      Tempf=abs(f(Solu(i)))
	  Tempx=Solu(i)
  Endif
 EndDo

  write(*,*)  "MAX |f(x)|=",tempf, "       ,Where x=", tempx

	 
contains

function f(x) result(f_result)
 real:: f_result,x
 f_result=2*x*cos(2*x)-(x-2)**2
end function f

function Df(x) result(Df_result)
 real:: Df_result,x
 Df_result=2*cos(2*x)-4*x*sin(2*x)-(x-2)*2
end function Df

function Df1(x) result(Df1_result)
 type(interval):: Df1_result,x
 Df1_result=2*cos(2*x)-4*x*sin(2*x)-(x-2)*2
end function Df1

function DDf(x) result(DDf_result)
 type(interval):: DDf_result,x
 DDf_result=-8*sin(2*x)-8*x*cos(2*x)-2
end function DDf


end program problem_3c_with_interval