bisect.f

c      PROGRAM TO FIND A ROOT OF  f(x) = 0  USING BISECTION
         program bisect
       implicit none
       real a,b,tol
       real ff
       real root,cc,y
       integer iter,itstart,itend

c      THE INTERVAL CONTAINING ROOT IS [a,b]
c      FIND ROOT TO WITHIN A TOLERANCE  tol
c      THE FUNCTION f(x) IS DEFINED IN A FUNCTION SUBROUTINE GIVING ff
       
c      DO A MAXIMUM OF itend ITERATIONS
       itstart=1
       itend=50

c      ENTER INTERVAL ENDPOINTS
       write(*,*) "Need interval [a,b] which contains root"
       write(*,*) "What is a?"
       read(*,*)a
       write(*,*) "What is b??"
       read(*,*)b
       write(*,*) "f(a) = ",ff(a)
       write(*,*) "f(b) = ",ff(b)
       write(*,*) "To what tolerance?"
       read(*,*)tol

c      CHECK THAT a IS NOT A ROOT
       if (abs(ff(a)) .lt. tol) then
         root=a
         go to 666
       endif
c      SIMILARLY FOR b
       if (abs(ff(b)) .lt. tol) then
         root=b
         go to 666
       endif
c      CHECK THAT f(a) AND f(b) HAVE OPPOSITE SIGNS
       if (ff(a)*ff(b) .ge. 0.0) then
         write(*,*) "f(a) and f(b) have same signs....try again"
         go to 999
       endif

c      INTERCHANGE a AND b TO ENSURE  a < b
       if (b .lt. a) then
         cc=a
         a=b
         b=cc
       endif

C      MAIN ITERATION LOOP
       do iter=itstart,itend

         write(*,*)"Interval is [",a,b,"]"

c        GET THE MIDPOINT
         y=(a+b)/2.0
c        IS y A GOOD ESTIMATE OF THE ROOT?
         if (abs(ff(y)) .lt. tol) then
           root = y
           go to 666
         endif
c        IF NOT, DETERMINE WHICH ENDPOINT TO KEEP
         if (ff(y)*ff(a) .lt. 0) then
           b=y
          else
           a=y
         endif
c        HAVE WE SHRUNK THE INTERVAL TOO SMALL?
         if (b-a .lt. tol) then
           write(*,*)"Interval is too small"
           write(*,*)"Interval is [",a,b,"]"
           go to 999
         endif

       enddo

666    write(*,*) "The root is",root," and f(y) = ",ff(root)

999    stop
       end

       real function ff(x)
c      THIS FUNCTION SUBROUTINE COMPUTES THE VALUE OF f AT A POINT x
       implicit none
       real x,ff

       ff=x*x-2.0

       return
       end