************************************************************************
*   FUNCTION:  Driver for using zero finding routines.                 *
*   AUTHORS:   Lawrence Shampine, Richard Allen, Steven Pruess  for    *
*              the text  Fundamentals of Numerical Computing           *
*   DATE:      June 12, 1995                                           *
************************************************************************
C     Sample program to illustrate ZERO.
C
      INTEGER          FLAG,NFEVAL
      DOUBLE PRECISION B,C,ABSERR,RELERR,F,RESIDL
C
C     Define labelled COMMON for the global variable NFEVAL.
C
      COMMON          /COST/ NFEVAL
      EXTERNAL         F
      EXTERNAL         F1 
C
C     Assign bracket and error tolerances.
C
      B = 1.D0
c     C = 0.5D0
      C = 0.38268343236509d0
      ABSERR = 1.D-8
      RELERR = 1.D-6
C
C     Initialize the count of the number of calls to F from ZERO.
C
      NFEVAL = 0
C
C     Use ZERO to find a root.
C
      CALL ZERO(F1,B,C,ABSERR,RELERR,RESIDL,FLAG)
C
C     Check FLAG and print results.
C
      IF (FLAG .EQ. -2) THEN
         PRINT *,' Initial interval [B,C] does not bracket a root.'
      ELSEIF (FLAG .EQ. -1) THEN
         PRINT *,' Error in specifying tolerances.'
      ELSEIF (FLAG .EQ. 0) THEN
         PRINT *,' Computed a root B = ',B
         PRINT *,' The residual F(B) = ',RESIDL
         PRINT *,NFEVAL,' evaluations of F were required.'
      ELSEIF (FLAG .EQ. 1) THEN
         PRINT *,' Too much work: NFEVAL = ',NFEVAL
         PRINT *,' There is a root in [B,C] with'
         PRINT *,' B = ',B,' C = ',C
         PRINT *,' The residual F(B)  = ',RESIDL
      ELSEIF (FLAG .EQ. 2) THEN
         PRINT *,' Computed a pole B = ',B
         PRINT *,' The residual F(B)  = ',RESIDL
      ENDIF
      STOP
      END
C-----------------------------------------------------------------------
      DOUBLE PRECISION FUNCTION F(X)
      DOUBLE PRECISION X
      INTEGER          NFEVAL
      COMMON          /COST/ NFEVAL
C
      F = EXP(-X) - 2.D0*X
      NFEVAL = NFEVAL+1
      RETURN
      END
C-----------------------------------------------------------------------
      DOUBLE PRECISION FUNCTION F1(X)
      IMPLICIT NONE
      DOUBLE PRECISION X
      DOUBLE PRECISION A(0:100)
      DOUBLE PRECISION P(0:1)
      INTEGER          NFEVAL
      COMMON          /COST/ NFEVAL
C
      A(0) = 128.d0
      A(1) =   0.d0
      A(2) =-256.d0
      A(3) =   0.d0
      A(4) = 160.d0
      A(5) =   0.d0
      A(6) = -32.d0
      A(7) =   0.d0
      A(8) =   1.d0
      CALL POLY(8,0,A,X,P)
      F1 = P(0)
      NFEVAL = NFEVAL+1
      RETURN
      END
C-----------------------------------------------------------------------
      DOUBLE PRECISION FUNCTION F2(X)
      IMPLICIT NONE
      DOUBLE PRECISION X
      DOUBLE PRECISION A(0:100)
      DOUBLE PRECISION P(0:1)
      INTEGER          NFEVAL
      COMMON          /COST/ NFEVAL
C
      A(0) = 256.d0
      A(1) =   0.d0
      A(2) =-576.d0
      A(3) =   0.d0
      A(4) = 432.d0
      A(5) =   0.d0
      A(6) =-120.d0
      A(7) =   0.d0
      A(8) =   9.d0
      A(9) =   0.d0
      CALL POLY(9,0,A,X,P)
      F1 = P(0)
      NFEVAL = NFEVAL+1
      RETURN
      END
C----------------------------------------------------------------------
        subroutine poly(n,k,A,X,P)
C computes the value of a polynomial of degree n at x by Horner's rule.
C if the derivative at x is also desired (k=1), then perform that
C evaluation simultaneously and return the result in P(1), 
C at cost comparable to that for computing P(0)
	implicit none
	double precision A(0:100)
        DOUBLE PRECISION P(0:1)
        double precision x
	INTEGER n,k,i,j
        DO J = 0,K
          P(J) = A(0)
        END DO
        DO I = 1,N-k
            P(0) = A(i) + x*P(0)
            P(1) = P(0) + x*P(1)
        END DO
        DO I = N-k+1,N
            P(0) = A(N) + x*P(0)
        END DO
        write(*,*) 'values of poly and ',k,' derivatives at ',x 
        write(*,*),(P(i), i=0,k)
        return
        end
C----------------------------------------------------------------------
