      SUBROUTINE INTRP(X,Y,XOUT,YOUT,YPOUT,NEQN,KOLD,PHI,PSI)
C
C     SANDIA MATHEMATICAL PROGRAM LIBRARY
C     APPLIED MATHEMATICS DIVISION 2642
C     SANDIA LABORATORIES
C     ALBUQUERQUE, NEW MEXICO  87115
C     CONTROL DATA 6600 VERSION 6.1  JANUARY 1976
C   * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C   *                 ISSUED BY SANDIA LABORATORIES,
C   *                   A PRIME CONTRACTOR TO THE
C   *   UNITED STATES ENERGY RESEARCH AND DEVELOPMENT ADMINISTRATION
C   * * * * * * * * * * * * * *  NOTICE   * * * * * * * * * * * * * * *
C   * THIS REPORT WAS PREPARED AS AN ACCOUNT OF WORK SPONSORED BY THE
C   * UNITED STATES GOVERNMENT.  NEITHER THE UNITED STATES NOR THE
C   * UNITED STATES ENERGY RESEARCH AND DEVELOPMENT ADMINISTRATION,
C   * NOR ANY OF THEIR EMPLOYEES, NOR ANY OF THEIR CONTRACTORS,
C   * SUBCONTRACTORS, OR THEIR EMPLOYEES, MAKES ANY WARRANTY, EXPRESS
C   * OR IMPLIED, OR ASSUMES ANY LEGAL LIABILITY OR RESPONSIBILITY
C   * FOR THE ACCURACY, COMPLETENESS OR USEFULNESS OF ANY INFORMATION,
C   * APPARATUS, PRODUCT OR PROCESS DISCLOSED, OR REPRESENTS THAT ITS
C   * USE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS.
C   * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C   * THE PRIMARY DOCUMENT FOR THE LIBRARY OF WHICH THIS ROUTINE IS
C   * IS A PART IS SAND75-0545, PRINTED FEB 1976.
C   * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C   WRITTEN BY L. F. SHAMPINE AND M. K. GORDON
C
C   ABSTRACT
C
C
C   THE METHODS IN SUBROUTINE  STEP  APPROXIMATE THE SOLUTION NEAR  X
C   BY A POLYNOMIAL.  SUBROUTINE  INTRP  APPROXIMATES THE SOLUTION AT
C   XOUT  BY EVALUATING THE POLYNOMIAL THERE.  INFORMATION DEFINING THIS
C   POLYNOMIAL IS PASSED FROM  STEP  SO  INTRP  CANNOT BE USED ALONE.
C
C   THIS CODE IS COMPLETELY EXPLAINED AND DOCUMENTED IN THE TEXT,
C   COMPUTER SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS0  THE INITIAL
C   VALUE PROBLEM  BY L. F. SHAMPINE AND M. K. GORDON.
C
C   INPUT TO INTRP --
C
C   THE USER PROVIDES STORAGE IN THE CALLING PROGRAM FOR THE ARRAYS IN
C   THE CALL LIST
C      DIMENSION Y(NEQN),YOUT(NEQN),YPOUT(NEQN),PHI(NEQN,16),PSI(12)
C   AND DEFINES
C      XOUT -- POINT AT WHICH SOLUTION IS DESIRED.
C   THE REMAINING PARAMETERS ARE DEFINED IN  STEP  AND PASSED TO  INTRP
C   FROM THAT SUBROUTINE
C
C   OUTPUT FROM  INTRP --
C
C      YOUT(*) -- SOLUTION AT  XOUT
C      YPOUT(*) -- DERIVATIVE OF SOLUTION AT  XOUT
C   THE REMAINING PARAMETERS ARE RETURNED UNALTERED FROM THEIR INPUT
C   VALUES.  INTEGRATION WITH  STEP  MAY BE CONTINUED.
C
      DIMENSION Y(NEQN),YOUT(NEQN),YPOUT(NEQN),PHI(NEQN,16),PSI(12)
      DIMENSION G(13),W(13),RHO(13)
      DATA G(1)/1.0/,RHO(1)/1.0/
C
      HI = XOUT - X
      KI = KOLD + 1
      KIP1 = KI + 1
C
C   INITIALIZE W(*) FOR COMPUTING G(*)
C
      DO 5 I = 1,KI
        TEMP1 = I
 5      W(I) = 1.0/TEMP1
      TERM = 0.0
C
C   COMPUTE G(*)
C
      DO 15 J = 2,KI
        JM1 = J - 1
        PSIJM1 = PSI(JM1)
        GAMMA = (HI + TERM)/PSIJM1
        ETA = HI/PSIJM1
        LIMIT1 = KIP1 - J
        DO 10 I = 1,LIMIT1
 10       W(I) = GAMMA*W(I) - ETA*W(I+1)
        G(J) = W(1)
        RHO(J) = GAMMA*RHO(JM1)
 15     TERM = PSIJM1
C
C   INTERPOLATE
C
      DO 20 L = 1,NEQN
        YPOUT(L) = 0.0
 20     YOUT(L) = 0.0
      DO 30 J = 1,KI
        I = KIP1 - J
        TEMP2 = G(I)
        TEMP3 = RHO(I)
        DO 25 L = 1,NEQN
          YOUT(L) = YOUT(L) + TEMP2*PHI(L,I)
 25       YPOUT(L) = YPOUT(L) + TEMP3*PHI(L,I)
 30     CONTINUE
      DO 35 L = 1,NEQN
 35     YOUT(L) = Y(L) + HI*YOUT(L)
      RETURN
      END
