	PROGRAM D9R6

C	Driver for routine RTSEC

	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(*,20)

20	FORMAT(/1X,'Roots of F(x):')

	WRITE(*,30)

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

	DO 11 I=1,NB

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

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

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

40	FORMAT(1X,'Root',I2,2X,F12.6,E16.4)

11	CONTINUE

	END

C

	FUNCTION FNCT(X)

	FNCT=exp(x)-3.*x

	RETURN

	END

      FUNCTION RTSEC(FUNC,X1,X2,XACC)

      PARAMETER (MAXIT=30)

      FL=FUNC(X1)

      F=FUNC(X2)

      IF(ABS(FL).LT.ABS(F))THEN

	RTSEC=X1

	XL=X2

	SWAP=FL

	FL=F

	F=SWAP

      ELSE

	XL=X1

	RTSEC=X2

      ENDIF

      DO 11 J=1,MAXIT

	DX=(XL-RTSEC)*F/(F-FL)

	XL=RTSEC

	FL=F

	RTSEC=RTSEC+DX

	F=FUNC(RTSEC)

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

11    CONTINUE

      PAUSE 'RTSEC exceed maximum iterations'

      END

      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

