        PROGRAM D9R4

C       Driver for routine RTBIS

        EXTERNAL FNCT

        PARAMETER(N=100,NBMAX=20)

        DIMENSION XB1(NBMAX),XB2(NBMAX)

        WRITE(*,*)'X1=?,X2=?,NB=?'

        READ(*,*)X1,X2,NB

        CALL ZBRAK(FNCT,X1,X2,N,XB1,XB2,NB)

        WRITE(*,'(/1X,A)') 'Roots of FNCT:'

        WRITE(*,20)

20      FORMAT(/1X,T19,'x',T31,'F(x)'/)

        DO 11 I=1,NB

                XACC=(1.0E-6)*(XB1(I)+XB2(I))/2.0

                ROOT=RTBIS(FNCT,XB1(I),XB2(I),XACC)

                WRITE(*,30) I,ROOT,FNCT(ROOT)

30      FORMAT(1X,'root',I2,2X,F12.6,E16.4)

11      CONTINUE

        END

C

        FUNCTION FNCT(X)

        FNCT=EXP(X)-3.*X

        RETURN

        END

C

      FUNCTION RTBIS(FUNC,X1,X2,XACC)

      PARAMETER (JMAX=40)

      FMID=FUNC(X2)

      F=FUNC(X1)

      IF(F*FMID.GE.0.) PAUSE 'Root must be bracketed for bisection.'

      IF(F.LT.0.)THEN

        RTBIS=X1

        DX=X2-X1

      ELSE

        RTBIS=X2

        DX=X1-X2

      ENDIF

      DO 11 J=1,JMAX

        DX=DX*.5

        XMID=RTBIS+DX

        FMID=FUNC(XMID)

        IF(FMID.LE.0.)RTBIS=XMID

        IF(ABS(DX).LT.XACC .OR. FMID.EQ.0.) RETURN

11    CONTINUE

      PAUSE 'too many bisections'

      END

C

      SUBROUTINE ZBRAK(FX,X1,X2,N,XB1,XB2,NB)

      PARAMETER(NBMAX=20)

      DIMENSION XB1(NBMAX),XB2(NBMAX)

      NBB=NB

      NB=0

      X=X1

      DX=(X2-X1)/N

      FP=FX(X)

      DO 11 I=1,N

        X=X+DX

        FC=FX(X)

        IF(FC*FP.LT.0.) THEN

          NB=NB+1

          XB1(NB)=X-DX

          XB2(NB)=X

        ENDIF

        FP=FC

        IF(NBB.EQ.NB)RETURN

11    CONTINUE

      RETURN

      END

