differentiate.f

c      PROGRAM TO FIND THE DERIVATIVE OF  f(x) NUMERICALLY
c      AND COMPARE WITH THE ACTUAL DERIVATIVE
        program differentiate
       implicit none
       real*8 xeval,delta,xclose
       real*8 feval,fclose,ff,deriv,fprime,ffprime
       real*8 diff
       integer POW
       real*8 pi
c      COMMON BLOCK consts  CONTAINS CONSTANTS TO BE SHARED
       common/consts/pi
c       parameter (pi= 3.1415926535897931d0)


       pi= 3.1415926535897931d0

c      EVALUATE THE DERIVATIVE AT xeval
c      A NEARBY POINT IS xclose=xeval + delta
c      WHERE delta = 10^POW
c      feval,fclose ARE f AT THESE POINTS
c      deriv IS THE APPROXIMATE DERIVATIVE
c      fprime IS THE ACTUAL DERIVATIVE


c      ENTER POINT OF EVALUATION
       write(*,*) "At what point do you wish to evaluate
     1     the derivative?"
       read(*,*)xeval

       write(*,*)"POW  x          approx f'         exact f'
     1   difference"

c      EXPONENTIATE TO POWER POW
       do POW=-1,-16,-1

         delta=10.0d0**(POW)
         xclose=xeval+delta
         feval=ff(xeval)
         fclose=ff(xclose)
         deriv=(fclose-feval)/(xclose-xeval)
         fprime=ffprime(xeval)
         diff=fprime-deriv

         write(*,*)POW,xeval,deriv,fprime,diff

       enddo


       stop
       end

c****6****************************************************************72
       real*8 function ff(x)
c      THIS FUNCTION SUBROUTINE COMPUTES THE VALUE OF f AT A POINT x
       implicit none
       real*8 pi,x
       common/consts/pi
c       parameter (pi= 3.1415926535897931d0)

       ff = dsin(2.d0*pi*x)

       return
       end

c****6****************************************************************72
       real*8 function ffprime(x)
c      THIS FUNCTION SUBROUTINE COMPUTES THE VALUE OF f AT A POINT x
       implicit none
       real*8 pi,x
       common/consts/pi
c       parameter (pi= 3.1415926535897931d0)

       ffprime=2.d0*pi*dcos(2.d0*pi*x)

       return
       end
c****6****************************************************************72