      PROGRAM NUMDERIV
C     NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994
C     To accompany the text:
C     NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
C     Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
C     This free software is complements of the author.
C
C     Algorithm 6.1 (Differentiation Using Limits).
C     Section 6.1, Approximating the Derivative, Page 326
C
      PARAMETER(Tol=1E-7,Max=100)
      INTEGER M
      REAL DX,H,X
      DIMENSION D(0:Max),DX(0:Max),E(0:Max)
      CHARACTER ANS*1
      EXTERNAL F
10    CALL INPUTS(X)
      CALL DLIMIT(F,D,DX,E,X,H,M,Tol)
      CALL RESULT(D,DX,E,X,H,M)
      WRITE(9,*)' '
      WRITE(9,*) 'WANT TO COMPUTE ANOTHER VALUE OF F`(X) ? '
      READ(9,'(A)') ANS
      IF (ANS.EQ.'Y' .OR. ANS.EQ.'y') GOTO 10
      STOP
      END

      REAL FUNCTION F(X)
      REAL X
      F=COS(X)
      RETURN
      END

      SUBROUTINE PRINTFUN
      WRITE(9,*)'    F(X)  =  COS(X)'
      RETURN
      END

      SUBROUTINE DLIMIT(F,D,DX,E,X,H,M,Tol)
      PARAMETER(Max=100)
      INTEGER M,N
      REAL D,DX,E,H,R,X
      DIMENSION D(0:Max),DX(0:Max),E(0:Max),R(0:Max)
      CHARACTER ANS*1
      EXTERNAL F
      Small=1E-9
      H = 1
      DX(0) = H
      D(0) = 0.5*(F(X+H)-F(X-H))/H
      DO N=1,2
        H = H/2
        DX(N) = H
        D(N) = 0.5*(F(X+H)-F(X-H))/H
        E(N) = ABS(D(N)-D(N-1))
        R(N) = 2*E(N)/(ABS(D(N))+ABS(D(N-1))+Small)
      ENDDO
      N = 1
      WHILE ((E(N).GT.E(N+1) .OR. R(N).GT.Tol).AND.(N.LT.Max))
        H = H/2
        DX(N+2) = H
        D(N+2) = 0.5*(F(X+H)-F(X-H))/H
        E(N+2) = ABS(D(N+2)-D(N+1))
        R(N+2) = 2*E(N+2)/(ABS(D(N+2))+ABS(D(N+1))+Small)
        N = N+1
      REPEAT
      M=N
      RETURN
      END

      SUBROUTINE XLIMIT(F,D,DX,E,X,H,M,Tol)
C     This subroutine uses labeled DO loop(s).
      PARAMETER(Max=100)
      INTEGER M,N
      REAL D,DX,E,H,R,X
      DIMENSION D(0:Max),DX(0:Max),E(0:Max),R(0:Max)
      CHARACTER ANS*1
      EXTERNAL F
      Small=1E-9
      H = 1
      DX(0) = H
      D(0) = 0.5*(F(X+H)-F(X-H))/H
      DO 10 N=1,2
        H = H/2
        DX(N) = H
        D(N) = 0.5*(F(X+H)-F(X-H))/H
        E(N) = ABS(D(N)-D(N-1))
        R(N) = 2*E(N)/(ABS(D(N))+ABS(D(N-1))+Small)
10    CONTINUE
      N = 1
20    IF ((E(N).GT.E(N+1) .OR. R(N).GT.Tol).AND.(N.LT.Max)) THEN
        H = H/2
        DX(N+2) = H
        D(N+2) = 0.5*(F(X+H)-F(X-H))/H
        E(N+2) = ABS(D(N+2)-D(N+1))
        R(N+2) = 2*E(N+2)/(ABS(D(N+2))+ABS(D(N+1))+Small)
        N = N+1
        GOTO 20
      ENDIF
      M=N
      RETURN
      END

      SUBROUTINE INPUTS(X)
      INTEGER I
      REAL X
      DO 10 I=1,16
        WRITE(9,*)' '
10    CONTINUE
      WRITE(9,*)'THE LIMIT OF THE SYMMETRIC DIFFERENCE QUOTIENTS'
      WRITE(9,*)' '
      WRITE(9,*)'IS USED TO FIND AN APPROXIMATION TO  F`(X):'
      WRITE(9,*)' '
      WRITE(9,*)' '
      WRITE(9,*)'   F`(X)  =  LIM  [F(X+H) - F(X-H)]/2H'
      WRITE(9,*)'            H -> 0 '
      WRITE(9,*)' '
      WRITE(9,*)' '
      WRITE(9,*)'THE GIVEN FUNCTION IS:'
      WRITE(9,*)' '
      WRITE(9,*)' '
      CALL PRINTFUN
      WRITE(9,*)' '
      WRITE(9,*)' '
      WRITE(9,*)'ENTER THE VALUE  X = '
      READ(9,*) X
      RETURN
      END

      SUBROUTINE RESULT(D,DX,E,X,H,M)
      PARAMETER(Max=100)
      INTEGER I,M
      REAL B,D,DL,DR,DX,E,H,X
      DIMENSION D(0:Max),DX(0:Max),E(0:Max)
      EXTERNAL F
      WRITE(9,*)' '
      WRITE(9,*)'          H
     &                      ( F(X+H) - F(X-H) )/2H'
      DO 10 N=0,M-1
      WRITE(9,1000) DX(N),D(N)
10    CONTINUE
      WRITE(9,*)' '
      PAUSE
      WRITE(9,*)' '
      WRITE(9,*)'THE LIMIT OF THE SYMMETRIC DIFFERENCE QUOTIENTS WAS
     + FOUND'
      WRITE(9,*)' '
      WRITE(9,*)'   F`(X)  =  LIM  [F(X+H) - F(X-H)]/2H'
      WRITE(9,*)'            H -> 0 '
      WRITE(9,*)' '
      CALL PRINTFUN
      WRITE(9,*)' '
      WRITE(9,*)'   AFTER ',M,' ITERATIONS, H =',H
      WRITE(9,*)' '
      WRITE(9,*)'AN APPROXIMATION FOR THE VALUE OF THE DERIVATIVE IS:'
      WRITE(9,*)' '
      WRITE(9,*)'F`(',X,'  )  =',D(M-1)
      WRITE(9,*)' '
      WRITE(9,*)'        THE ACCURACY IS  +-',E(M-1)
      DL=(F(X)-F(X-H))/H
      DR=(F(X+H)-F(X))/H
      B=1+ABS(DL)+ABS(DR)
      IF (ABS(DR-DL)>0.001*B) THEN
        WRITE(9,*)' '
        WRITE(9,*)'THE LEFT-HAND AND RIGHT-HAND DERIVATIVES ARE NOT
     + EQUAL.'
        WRITE(9,*)' '
        WRITE(9,*)'LEFT-HAND   [F(X)-F(X-H)]/H =',DL
        WRITE(9,*)' '
        WRITE(9,*)'RIGHT-HAND  [F(X+H)-F(X)]/H =',DR
        WRITE(9,*)' '
        WRITE(9,*)'PERHAPS F`(X) DOES NOT EXIST.'
      ENDIF
      PAUSE
      RETURN
1000  FORMAT(F20.10,5X,F20.10)
      END

      PROGRAM DIFFEXTR
C     NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994
C     To accompany the text:
C     NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
C     Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
C     This free software is complements of the author.
C
C     Algorithm 6.2 (Differentiation Using Extrapolation).
C     Section 6.1, Approximating the Derivative, Page 327
C
      PARAMETER(Tol=1E-7,Delta=1E-7,MaxN=15)
      INTEGER J,K,N
      REAL D,DX,Error,H,RelErr,X
      DIMENSION D(0:MaxN,0:MaxN),DX(0:MaxN)
      CHARACTER ANS*1
      EXTERNAL F
10    CALL INPUTS(X)
      CALL EXTRAP(F,D,DX,X,H,Error,N,Tol,Delta)
      CALL RESULT(F,D,DX,X,H,Error,N)
      WRITE(9,*)' '
      WRITE(9,*)'WANT TO EVALUATE F`(X) AT ANOTHER POINT ? '
      READ(9,'(A)') ANS
      IF (ANS.EQ.'Y' .OR. ANS.EQ.'y') GOTO 10
      STOP
      END

      REAL FUNCTION F(X)
      REAL X
      F=SIN(X)
      RETURN
      END

      SUBROUTINE PRINTFUN
      WRITE(9,*)'   F(X)  =  SIN(X)'
      RETURN
      END

      SUBROUTINE EXTRAP(F,D,DX,X,H,Error,N,Tol,Delta)
      PARAMETER(MaxN=15)
      INTEGER J,K,N
      REAL D,DX,Error,H,RelErr,Small,X
      DIMENSION D(0:MaxN,0:MaxN),DX(0:MaxN)
      EXTERNAL F
      Small=1E-7
      H=1
      DX(0) = H
      J=1
      Error=1
      RelErr=1
      D(0,0)=0.5*(F(X+H)-F(X-H))/H
      WHILE ( RelErr.GT.Tol .AND. Error.GT.Delta .AND. J.LT.16)
        H=H/2
        DX(J) = H
        D(J,0)=0.5*(F(X+H)-F(X-H))/H
        DO K=1,J
          D(J,K)=D(J,K-1)+(D(J,K-1)-D(J-1,K-1))/(4**K-1)
        ENDDO
        Error=ABS(D(J,J)-D(J-1,J-1))
        RelErr=2*Error/(ABS(D(J,J))+ABS(D(J-1,J-1))+Small)
        N=J
        J=J+1
      REPEAT
      RETURN
      END

      SUBROUTINE XEXTRAP(F,D,DX,X,H,Error,N,Tol,Delta)
C     This subroutine uses simulated WHILE loop(s).
      PARAMETER(MaxN=15)
      INTEGER J,K,N
      REAL D,DX,Error,H,RelErr,Small,X
      DIMENSION D(0:MaxN,0:MaxN),DX(0:MaxN)
      EXTERNAL F
      Small=1E-7
      H=1
      DX(0) = H
      J=1
      Error=1
      RelErr=1
      D(0,0)=0.5*(F(X+H)-F(X-H))/H
10    IF ( RelErr.GT.Tol .AND. Error.GT.Delta .AND. J.LT.16) THEN
      H=H/2
      DX(J) = H
      D(J,0)=0.5*(F(X+H)-F(X-H))/H
      DO 20 K=1,J
        D(J,K)=D(J,K-1)+(D(J,K-1)-D(J-1,K-1))/(4**K-1)
20    CONTINUE
        Error=ABS(D(J,J)-D(J-1,J-1))
        RelErr=2*Error/(ABS(D(J,J))+ABS(D(J-1,J-1))+Small)
        N=J
        J=J+1
        GOTO 10
      ENDIF
      RETURN
      END

      SUBROUTINE INPUTS(X)
      INTEGER I
      REAL X
      DO 10 I=1,18
      WRITE(9,*)' '
10    CONTINUE
      WRITE(9,*)'RICHARDSON`S EXTRAPOLATION APPLIED TO THE SEQUENCE OF
     + DIFFERENCE QUOTIENTS'
      WRITE(9,*)' '
      WRITE(9,*)'                                                -K  '
      WRITE(9,*)'        { [F(X+H ) - F(X-H )]/2H : WHERE  H  = 2   }'
      WRITE(9,*)'                K         K     K          K'
      WRITE(9,*)' '
      WRITE(9,*)'IS USED TO FIND AN APPROXIMATION TO THE DERIVATIVE OF
     + THE FUNCTION:'
      WRITE(9,*)' '
      WRITE(9,*)' '
      CALL PRINTFUN
      WRITE(9,*)' '
      WRITE(9,*)' '
      WRITE(9,*)'ENTER THE VALUE  X = '
      READ(9,*) X
      RETURN
      END

      SUBROUTINE RESULT(F,D,DX,X,H,Error,N)
      PARAMETER(MaxN=15)
      INTEGER I,N
      REAL B,D,DL,DR,DX,Error,H,H0,X
      DIMENSION D(0:MaxN,0:MaxN),DX(0:MaxN)
      EXTERNAL F
      WRITE(9,*)' '
      WRITE(9,*)'          H                        D(N,N)'
      DO 10 I=0,N
      WRITE(9,1000) DX(I),D(I,I)
10    CONTINUE
      WRITE(9,*)' '
      WRITE(9,*)'RICHARDSON`S EXTRAPOLATION WAS USED ON THE SEQUENCE OF
     + DIFFERENCE QUOTIENTS'
      WRITE(9,*)'                                                    -K'
      WRITE(9,*)'          { [F(X+H ) - F(X-H )]/2H  : WHERE  H  = 2
     +}'
      WRITE(9,*)'                  K         K     K           K'
      WRITE(9,*)'TO FIND AN APPROXIMATION TO THE DERIVATIVE OF THE
     + FUNCTION:'
      WRITE(9,*)' '
      CALL PRINTFUN
      WRITE(9,*)' '
      WRITE(9,*)'   AFTER ',N,' ITERATIONS, H =',H
      WRITE(9,*)' '
      WRITE(9,*)'AN APPROXIMATION FOR THE VALUE OF THE DERIVATIVE IS'
      WRITE(9,*)' '
      WRITE(9,*)'F`(',X,'  )  =',D(N,N)
      WRITE(9,*)' '
      WRITE(9,*)'        THE ACCURACY IS  +-',Error
      H0=0.002
      IF (H.LT.H0) THEN
        H0=H
      ENDIF
      DL=(F(X)-F(X-H0))/H0
      DR=(F(X+H0)-F(X))/H0
      B=1+ABS(DL)+ABS(DR)
      IF (ABS(DR-DL).GT.0.02*B) THEN
        WRITE(9,*)' '
        WRITE(9,*)'THE LEFT-HAND AND RIGHT-HAND DERIVATIVES ARE NOT
     + EQUAL.'
        WRITE(9,*)' '
        WRITE(9,*)'LEFT-HAND   [F(X)-F(X-H0)]/H0 =',DL
        WRITE(9,*)' '
        WRITE(9,*)'RIGHT-HAND  [F(X+H0)-F(X)]/H0 =',DR
        WRITE(9,*)' '
        WRITE(9,*)'PERHAPS F`(X) DOES NOT EXIST.'
      ENDIF
      PAUSE
      RETURN
1000  FORMAT(F20.10,5X,F20.10)
      END

      PROGRAM NEWPOLYD
C     NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994
C     To accompany the text:
C     NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
C     Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
C     This free software is complements of the author.
C
C     Algorithm 6.3 (Differentiation Based on N+1 Nodes).
C     Section 6.2, Numerical Differentiation Formulas, Page 342
C
      PARAMETER(MaxN=15)
      INTEGER N
      REAL A,Df,X,Y
      CHARACTER ANS*1
      DIMENSION A(0:MaxN),X(0:MaxN),Y(0:MaxN)
10    CALL GETDATA(X,Y,N)
      CALL DIFPOLY(A,X,Y,N,Df)
      CALL RESULTS(A,X,Y,N,Df)
      WRITE(9,*)' '
      WRITE(9,*)'WANT TO USE A DIFFERENT SET OF DATA POINTS ? '
      READ(9,'(A)') ANS
      IF (ANS.EQ.'Y' .OR. ANS.EQ.'y') GOTO 10
      END

      SUBROUTINE DIFPOLY(A,X,Y,N,Df)
      PARAMETER(MaxN=15)
      INTEGER J,K,N
      REAL A,Df,Prod,T,X,Y
      DIMENSION A(0:MaxN),X(0:MaxN),Y(0:MaxN)
      DO K=0,N
        A(K)=Y(K)
      ENDDO
      DO J=1,N
        DO K=N,J,-1
          A(K)=(A(K)-A(K-1))/(X(K)-X(K-J))
        ENDDO
      ENDDO
      T=X(0)
      Df=A(1)
      Prod=1
      DO K=2,N
        Prod=Prod*(T-X(K-1))
        Df=Df+Prod*A(K)
      ENDDO
      RETURN
      END

      SUBROUTINE XDIFPOLY(A,X,Y,N,Df)
C     This subroutine uses labeled DO loop(s).
      PARAMETER(MaxN=15)
      INTEGER J,K,N
      REAL A,Df,Prod,T,X,Y
      DIMENSION A(0:MaxN),X(0:MaxN),Y(0:MaxN)
      DO 10 K=0,N
        A(K)=Y(K)
10    CONTINUE
      DO 30 J=1,N
      DO 20 K=N,J,-1
        A(K)=(A(K)-A(K-1))/(X(K)-X(K-J))
20    CONTINUE
30    CONTINUE
      T=X(0)
      Df=A(1)
      Prod=1
      DO 40 K=2,N
        Prod=Prod*(T-X(K-1))
        Df=Df+Prod*A(K)
40    CONTINUE
      RETURN
      END

      SUBROUTINE PRINTPOL(A,X,Y,N)
      PARAMETER(MaxN=15)
      INTEGER K,N
      REAL A,X,Y
      DIMENSION A(0:MaxN),X(0:MaxN),Y(0:MaxN)
        WRITE(9,*)' '
      IF (N.EQ.1) THEN
        WRITE(9,*)'P(T) = A  + A [T-X ]'
        WRITE(9,*)'        0    1    0 '
      ELSEIF (N.EQ.2) THEN
        WRITE(9,*)'P(T) = A  + A [T-X ] + A [T-X ][T-X ]'
        WRITE(9,*)'        0    1    0     2    0     1 '
      ELSEIF (N.EQ.3) THEN
        WRITE(9,*)'P(T) = A  + A [T-X ] + A [T-X ][T-X ] + A [T-X ][T-X
     +][T-X ]'
        WRITE(9,*)'        0    1    0     2    0     1     3    0     1
     +     2 '
      ELSE
        WRITE(9,*)'P(T) = A  + A [T-X ] + A [T-X ][T-X ] + A [T-X ][T-X
     +][T-X ]'
        WRITE(9,*)'        0    1    0     2    0     1     3    0     1
     +     2'
        WRITE(9,*)'       +...+ A [T-X ][T-X ] ... [T-X ]'
        WRITE(9,*)'    ',N,'    0     1',(N-1)
      ENDIF
      WRITE(9,*)' '
      WRITE(9,*)'         THE COEFFICIENTS:  THE ABSCISSAS:       THE OR
     +DINATES:'
      WRITE(9,*)'K            A(K)                X(K)                Y(
     +K)'
      WRITE(9,*)' '
      DO 10 K=0,N
        WRITE(9,999) K,A(K),X(K),Y(K)
10    CONTINUE
999   FORMAT(1X,I2,3F20.7)
      RETURN
      END

      SUBROUTINE GETDATA(X,Y,N)
      PARAMETER(MaxN=15)
      INTEGER I,K,N
      REAL X,Y
      DIMENSION X(0:MaxN),Y(0:MaxN)
      DO 10 I=1,18
      WRITE(9,*)' '
10    CONTINUE
      WRITE(9,*)'     THE DERIVATIVE OF THE NEWTON POLYNOMIAL'
      WRITE(9,*)' '
      WRITE(9,*)'P(T) = A  + A [T-X ] + A [T-X ][T-X ] + A [T-X ][T-X ][
     +T-X ]'
      WRITE(9,*)'        0    1    0     2    0     1     3    0     1
     +   2'
      WRITE(9,*)'       +...+ A [T-X ][T-X ]...[T-X   ]'
      WRITE(9,*)'              N    0     1        N-1'
      WRITE(9,*)' '
      WRITE(9,*)'IS USED TO APPROXIMATE THE DERIVATIVE AT THE POINT (X ,
     +Y ).'
      WRITE(9,*)'                                                     0
     + 0'
      WRITE(9,*)' '
      WRITE(9,*)'P(T) IS BASED ON THE N+1 NODES (X ,Y ),(X ,Y ),...,(X ,
     +Y ).'
      WRITE(9,*)'                                 0  0    1  1        N
     + N  '
      WRITE(9,*)' '
      WRITE(9,*)'THE VALUES {X } DO NOT NEED TO BE INCREASING.'
      WRITE(9,*)'             K'
      WRITE(9,*)' '
      WRITE(9,*)'ENTER THE NUMBER  N = '
      READ(9,*) N
      WRITE(9,*)' '
      WRITE(9,*)'NOW ENTER THE ',N+1,' POINTS'
      WRITE(9,*)' '
      DO 20 K=0,N
        WRITE(9,*)'X(',K,') = '
        READ(9,*) X(K)
        WRITE(9,*)'Y(',K,') = '
        READ(9,*) Y(K)
        WRITE(9,*)' '
20    CONTINUE
      RETURN
      END

      SUBROUTINE RESULTS(A,X,Y,N,Df)
      PARAMETER(MaxN=15)
      INTEGER I,N
      REAL A,Df,X,Y
      DIMENSION A(0:MaxN),X(0:MaxN),Y(0:MaxN)
      DO 10 I=1,18
        WRITE(9,*)' '
10    CONTINUE
      WRITE(9,*)'    THE NEWTON POLYNOMIAL OF DEGREE ',N,' IS:'
      WRITE(9,*)' '
      CALL PRINTPOL(A,X,Y,N)
      WRITE(9,*)' '
      WRITE(9,*)'THE APPROXIMATION FOR THE DERIVATIVE AT X  =',X(0)
      WRITE(9,*)'                                         0'
      WRITE(9,*)'USING THE NEWTON POLYNOMIAL IS:'
      WRITE(9,*)' '
      WRITE(9,*)'P`(',X(0),'  ) = ',Df
      WRITE(9,*)' '
      RETURN
      END