1      SUBROUTINE INTRP(X,Y,XOUT,YOUT,YPOUT,NEQN,KOLD,PHI,PSI)
2C***BEGIN PROLOGUE  INTRP
3C***DATE WRITTEN   740101   (YYMMDD)
4C***REVISION DATE  820801   (YYMMDD)
5C***CATEGORY NO.  Z
6C***KEYWORDS  ADAMS METHOD,DEPAC,INITIAL VALUE PROBLEMS,ODE,
7C             ORDINARY DIFFERENTIAL EQUATIONS,PREDICTOR-CORRECTOR
8C***AUTHOR  SHAMPINE, L. F., (SNLA)
9C           GORDON, M. K., (SNLA)
10C***PURPOSE  Approximates the solution at XOUT by evaluating the
11C            polynomial computed in STEP2 at XOUT.  Must be used in
12C            conjunction with STEP2.
13C***DESCRIPTION
14C
15C   Written by L. F. Shampine and M. K. Gordon
16C
17C   Abstract
18C
19C
20C   The methods in subroutine  STEP2  approximate the solution near  X
21C   by a polynomial.  Subroutine  INTRP  approximates the solution at
22C   XOUT  by evaluating the polynomial there.  Information defining this
23C   polynomial is passed from  STEP2  so  INTRP  cannot be used alone.
24C
25C   This code is completely explained and documented in the text,
26C   "Computer Solution of Ordinary Differential Equations, the Initial
27C   Value Problem"  by L. F. Shampine and M. K. Gordon.
28C   Further details on use of this code are available in "Solving
29C   Ordinary Differential Equations with ODE, STEP, and INTRP",
30C   by L. F. Shampine and M. K. Gordon, SLA-73-1060.
31C
32C   Input to INTRP --
33C
34C   The user provides storage in the calling program for the arrays in
35C   the call list
36C      DIMENSION Y(NEQN),YOUT(NEQN),YPOUT(NEQN),PHI(NEQN,16),PSI(12)
37C   and defines
38C      XOUT -- point at which solution is desired.
39C   The remaining parameters are defined in  STEP2  and passed to
40C   INTRP  from that subroutine
41C
42C   Output from  INTRP --
43C
44C      YOUT(*) -- solution at  XOUT
45C      YPOUT(*) -- derivative of solution at  XOUT
46C   The remaining parameters are returned unaltered from their input
47C   values.  Integration with  STEP2  may be continued.
48C***REFERENCES  SHAMPINE, L. F., GORDON, M. K., *SOLVING
49C                 ORDINARY DIFFERENTIAL EQUATIONS WITH ODE, STEP,
50C                 AND INTRP*, SLA-73-1060, SANDIA LABORATORIES,
51C                 1973.
52C***ROUTINES CALLED  (NONE)
53C***END PROLOGUE  INTRP
54C
55C
56      DIMENSION Y(NEQN),YOUT(NEQN),YPOUT(NEQN),PHI(NEQN,16),PSI(12)
57      DIMENSION G(13),W(13),RHO(13)
58      DATA G(1)/1.0/,RHO(1)/1.0/
59C
60C***FIRST EXECUTABLE STATEMENT  INTRP
61      HI = XOUT - X
62      KI = KOLD + 1
63      KIP1 = KI + 1
64C
65C   INITIALIZE W(*) FOR COMPUTING G(*)
66C
67      DO 5 I = 1,KI
68        TEMP1 = I
69 5      W(I) = 1.0/TEMP1
70      TERM = 0.0
71C
72C   COMPUTE G(*)
73C
74      DO 15 J = 2,KI
75        JM1 = J - 1
76        PSIJM1 = PSI(JM1)
77        GAMMA = (HI + TERM)/PSIJM1
78        ETA = HI/PSIJM1
79        LIMIT1 = KIP1 - J
80        DO 10 I = 1,LIMIT1
81 10       W(I) = GAMMA*W(I) - ETA*W(I+1)
82        G(J) = W(1)
83        RHO(J) = GAMMA*RHO(JM1)
84 15     TERM = PSIJM1
85C
86C   INTERPOLATE
87C
88      DO 20 L = 1,NEQN
89        YPOUT(L) = 0.0
90 20     YOUT(L) = 0.0
91      DO 30 J = 1,KI
92        I = KIP1 - J
93        TEMP2 = G(I)
94        TEMP3 = RHO(I)
95        DO 25 L = 1,NEQN
96          YOUT(L) = YOUT(L) + TEMP2*PHI(L,I)
97 25       YPOUT(L) = YPOUT(L) + TEMP3*PHI(L,I)
98 30     CONTINUE
99      DO 35 L = 1,NEQN
100 35     YOUT(L) = Y(L) + HI*YOUT(L)
101      RETURN
102      END
103