      program bisektion
*---------------------------------------------------------------
*
*     Nullstellensuche mit Bisektion
*     
*---------------------------------------------------------------
      implicit none
*
      real*8     xlinks,xrechts,xmitte
      real*8     f
      real*8     flinks,frechts,fmitte
      real*8     Genauigkeit
*     
      character  Version*50
*
      Version= '# bisektion 01 Juergen Schnack  22. 05. 2009'
*---------------------------------------------------------------
*
*     Version
*
      write(6,'(a)') Version
*---------------------------------------------------------------
*     
*     Einlesen von xlinks, xrechts, Genauigkeit
*
      write(6,'(a)',advance="no") 'xlinks = '
      read(5,*) xlinks
      write(6,'(a)',advance="no") 'xrechts = '
      read(5,*) xrechts
      write(6,'(a)',advance="no") 'Genauigkeit = '
      read(5,*) Genauigkeit

*---------------------------------------------------------------
*     
*     Notwendige Voraussetzungen Prüfen
*
      flinks = f(xlinks)
      frechts = f(xrechts)
      IF((frechts*flinks).gt.0.d0) then
         write(6,'(a)') '# Ungültige Intervallgrenzen!'
         goto 9999
      endif
      xmitte = xlinks
      fmitte = flinks

*---------------------------------------------------------------
*     
*     Bisektion
*
      DO WHILE (DABS(fmitte).gt.Genauigkeit)
         xmitte = 0.5d0*(xrechts+xlinks)
         fmitte = f(xmitte)
         IF((fmitte*flinks).gt.0.d0) then
            xlinks = xmitte
            flinks = fmitte
         else
            xrechts = xmitte
            frechts = fmitte
         endif
         write(6,'(4f15.8)') xlinks,xrechts,fmitte,(xrechts-xlinks)
      END DO
      write(6,'(a,2f15.8)') '===> ',xmitte,fmitte
*---------------------------------------------------------------
*
*     Schluss
*
 9999 continue
      write(6,'(a)') '# Fertig!'

      stop
      end
*===============================================================
      real*8 function f(x)
      implicit none
      real*8   x

      f = (x+1.d0)*(x-2.d0)*(x-10.d0)

      return
      end

