ccccccccccccccccccccccccc     Program 2.4     cccccccccccccccccccccccccc
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                                                                      c
c Please Note:                                                         c
c                                                                      c
c (1) This computer program is part of the book, "An Introduction to   c
c     Computational Physics," written by Tao Pang and published and    c
c     copyrighted by Cambridge University Press in 1997.               c
c                                                                      c
c (2) No warranties, express or implied, are made for this program.    c
c                                                                      c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      PROGRAM DERIVATIVES
C
C Main program for derivatives of f(x) = sin(x).  F1: f';
C F2: f";  D1: error in f'; and D2: error in f".
C
      PARAMETER (N = 101)
      DIMENSION X(N),F(N),F1(N),D1(N),F2(N),D2(N)
C
      PI = 4.0*ATAN(1.0)
      H  = (PI/2./100
      DO      100  I = 1, N
        X(I) = H*(I-1)
        F(I) = SIN(X(I))
  100 CONTINUE
      CALL THREE(N,H,F,F1,F2)
      DO      300  I = 1, N
        D1(I) = F1(I)-COS(X(I))
        D2(I) = F2(I)+SIN(X(I))
        WRITE (6,999) X(I),F1(I),D1(I),F2(I),D2(I)
  300 CONTINUE
      STOP
  999 FORMAT (5F10.6)
      END
C
      SUBROUTINE THREE(N,H,FI,F1,F2)
C
C Subroutine for 1st and 2nd order derivatives with the three-
C point formulas. Extrapolations are made at the boundaries.
C FI: input f(x); H: interval; F1: f'; and F2: f".
C
      DIMENSION FI(N),F1(N),F2(N)
C
C f' and f" from three-point formulas
C
      DO      100  I = 2, N-1
        F1(I) = (FI(I+1) - FI(I-1))/(2.*H)
        F2(I) = (FI(I+1) - 2.*FI(I) + FI(I-1))/(H*H)
  100 CONTINUE
C
C Linear extrapolation for the boundary points
C
      F1(1) = 2.*F1(2) - F1(3)
      F1(N) = 2.*F1(N-1) - F1(N-2)
      F2(1) = 2.*F2(2) - F2(3)
      F2(N) = 2.*F2(N-1) - F2(N-2)
      RETURN
      END