C
C PAGE 300: NUMERICAL MATHEMATICS AND COMPUTING, CHENEY/KINCAID, 1985
C
C FILE: SCH.FOR
C
C INTERPOLATES TABLE USING SCHOENBERG'S PROCESS (SCH,ESCH,F)
C
      DIMENSION D(22),X(101),Y(101),Z(101)
      EXTERNAL F, ESCH
      DATA A,B,N,NP2,M/-5.,5.,20,22,101/
      CALL SCH(F,A,B,NP2,D) 
      PRINT 2,(D(J),J=1,NP2)
 2    FORMAT(1X,5E25.15)
      H=(B-A)/FLOAT(M-1)
      DO 3 J=1,M
      X(J)=A+H*FLOAT(J-1) 
      Y(J)=ESCH(A,B,NP2,D,X(J)) 
 3    Z(J)=Y(J)-F(X(J)) 
      DO 4 J=1,M
 4    PRINT 2,X(J),Y(J),Z(J)
      STOP
      END 
  
      FUNCTION F(X) 
      F=1.0/(X**2+1.0)
      RETURN
      END 
  
      SUBROUTINE SCH(F,A,B,NP2,D)     
      DIMENSION  D(NP2)     
      H = (B - A)/REAL(NP2-3) 
      DO 2 J = 2,NP2-1      
   2  D(J) = F(A + H*REAL(J-2))       
      D(1) = 2.0*D(2) - D(3)
      D(NP2) = 2.0*D(NP2-1) - D(NP2-2)
      RETURN      
      END 
  
      FUNCTION ESCH(A,B,NP2,D,X)      
      DIMENSION  D(NP2)     
      H = (B - A)/REAL(NP2-3) 
      C = A - 2.5*H 
      Y = (X - C)/H 
      J = INT(Y)  
      Y = Y - AINT(Y)       
      P = D(J+1)*Y + D(J)*(2.0 - Y)   
      Q = D(J)*(Y + 1.0) + D(J-1)*(1.0 - Y)     
      ESCH = (P*Y + Q*(1.0 - Y))*0.5  
      RETURN      
      END