c
c     Numerical Analysis:
c     The Mathematics of Scientific Computing
c     D.R. Kincaid & E.W. Cheney
c     Brooks/Cole Publ., 1990
c
c     Section 3.1
c
c     Example of Bisection Method
c
c
c     file: ex1s31.f 
c
      external f
      data a,b,epsi,M,delta/-4.0,-3.0,1.0e-5,25,1.0e-5/	
c
      print *
      print *,' Bisection Method to find the root of exp(x)=sin x'
      print *,' Section 3.1, Kincaid-Cheney'
      print *
c
      call bisect(f,a,b,epsi,delta,M)
      stop
      end 
c
      function f(x) 
      f =  exp(x) - sin(x)
      return
      end 
c
      subroutine bisect(f,a,b,epsi,delta,M)       
c
c     bisecting the interval
c
      external f
      fa = f(a)
      fb = f(b)
      error = b - a
      print *,'    a =',a,'     b =',b
      print *,' f(a) =',fa,'  f(b) =',fb
      print *
      if (sign(1.0,fa) .eq. sign(1.0,fb)) then
         print 5
         return
      endif
      print 3
c
      do 2 k=1,M
         error = error/2.0   
         c = a + error
         fc = f(c) 
         print 4,k,c,fc,error
         if ((abs(fc) .lt. epsi) .or. (abs(error) .lt. delta))  return
         if (sign(1.0,fa) .ne. sign(1.0,fc)) then 
            b = c 
            fb = fc 
         else
            a = c 
            fa = fc 
         endif 
 2    continue    
c
 3    format (2x,' step',6x,'c',13x,'f(c)',10x,'error')
 4    format (2x,i3,3(2x,e13.6))  
 5    format (2x,'function does not change sign over interval',2e13.6)
      return
      end