1*DECK SINTRP
2      SUBROUTINE SINTRP (X, Y, XOUT, YOUT, YPOUT, NEQN, KOLD, PHI, IVC,
3     +   IV, KGI, GI, ALPHA, OG, OW, OX, OY)
4C***BEGIN PROLOGUE  SINTRP
5C***PURPOSE  Approximate the solution at XOUT by evaluating the
6C            polynomial computed in STEPS at XOUT.  Must be used in
7C            conjunction with STEPS.
8C***LIBRARY   SLATEC (DEPAC)
9C***CATEGORY  I1A1B
10C***TYPE      SINGLE PRECISION (SINTRP-S, DINTP-D)
11C***KEYWORDS  ADAMS METHOD, DEPAC, INITIAL VALUE PROBLEMS, ODE,
12C             ORDINARY DIFFERENTIAL EQUATIONS, PREDICTOR-CORRECTOR,
13C             SMOOTH INTERPOLANT
14C***AUTHOR  Watts, H. A., (SNLA)
15C***DESCRIPTION
16C
17C   The methods in subroutine  STEPS  approximate the solution near  X
18C   by a polynomial.  Subroutine  SINTRP  approximates the solution at
19C   XOUT  by evaluating the polynomial there.  Information defining this
20C   polynomial is passed from  STEPS  so  SINTRP  cannot be used alone.
21C
22C   Subroutine STEPS is completely explained and documented in the text,
23C   "Computer Solution of Ordinary Differential Equations, the Initial
24C   Value Problem"  by L. F. Shampine and M. K. Gordon.
25C
26C   Input to SINTRP --
27C
28C   The user provides storage in the calling program for the arrays in
29C   the call list
30C      DIMENSION Y(NEQN),YOUT(NEQN),YPOUT(NEQN),PHI(NEQN,16),OY(NEQN)
31C                AND ALPHA(12),OG(13),OW(12),GI(11),IV(10)
32C   and defines
33C      XOUT -- point at which solution is desired.
34C   The remaining parameters are defined in  STEPS  and passed to
35C   SINTRP  from that subroutine
36C
37C   Output from  SINTRP --
38C
39C      YOUT(*) -- solution at  XOUT
40C      YPOUT(*) -- derivative of solution at  XOUT
41C   The remaining parameters are returned unaltered from their input
42C   values.  Integration with  STEPS  may be continued.
43C
44C***REFERENCES  H. A. Watts, A smoother interpolant for DE/STEP, INTRP
45C                 II, Report SAND84-0293, Sandia Laboratories, 1984.
46C***ROUTINES CALLED  (NONE)
47C***REVISION HISTORY  (YYMMDD)
48C   840201  DATE WRITTEN
49C   890831  Modified array declarations.  (WRB)
50C   890831  REVISION DATE from Version 3.2
51C   891214  Prologue converted to Version 4.0 format.  (BAB)
52C   920501  Reformatted the REFERENCES section.  (WRB)
53C***END PROLOGUE  SINTRP
54C
55      DIMENSION Y(*),YOUT(*),YPOUT(*),PHI(NEQN,16),OY(*)
56      DIMENSION G(13),C(13),W(13),OG(13),OW(12),ALPHA(12),GI(11),IV(10)
57C
58C***FIRST EXECUTABLE STATEMENT  SINTRP
59      KP1 = KOLD + 1
60      KP2 = KOLD + 2
61C
62      HI = XOUT - OX
63      H = X - OX
64      XI = HI/H
65      XIM1 = XI - 1.
66C
67C   INITIALIZE W(*) FOR COMPUTING G(*)
68C
69      XIQ = XI
70      DO 10 IQ = 1,KP1
71        XIQ = XI*XIQ
72        TEMP1 = IQ*(IQ+1)
73 10     W(IQ) = XIQ/TEMP1
74C
75C   COMPUTE THE DOUBLE INTEGRAL TERM GDI
76C
77      IF (KOLD .LE. KGI) GO TO 50
78      IF (IVC .GT. 0) GO TO 20
79      GDI = 1.0/TEMP1
80      M = 2
81      GO TO 30
82 20   IW = IV(IVC)
83      GDI = OW(IW)
84      M = KOLD - IW + 3
85 30   IF (M .GT. KOLD) GO TO 60
86      DO 40 I = M,KOLD
87 40     GDI = OW(KP2-I) - ALPHA(I)*GDI
88      GO TO 60
89 50   GDI = GI(KOLD)
90C
91C   COMPUTE G(*) AND C(*)
92C
93 60   G(1) = XI
94      G(2) = 0.5*XI*XI
95      C(1) = 1.0
96      C(2) = XI
97      IF (KOLD .LT. 2) GO TO 90
98      DO 80 I = 2,KOLD
99        ALP = ALPHA(I)
100        GAMMA = 1.0 + XIM1*ALP
101        L = KP2 - I
102        DO 70 JQ = 1,L
103 70       W(JQ) = GAMMA*W(JQ) - ALP*W(JQ+1)
104        G(I+1) = W(1)
105 80     C(I+1) = GAMMA*C(I)
106C
107C   DEFINE INTERPOLATION PARAMETERS
108C
109 90   SIGMA = (W(2) - XIM1*W(1))/GDI
110      RMU = XIM1*C(KP1)/GDI
111      HMU = RMU/H
112C
113C   INTERPOLATE FOR THE SOLUTION -- YOUT
114C   AND FOR THE DERIVATIVE OF THE SOLUTION -- YPOUT
115C
116      DO 100 L = 1,NEQN
117        YOUT(L) = 0.0
118 100    YPOUT(L) = 0.0
119      DO 120 J = 1,KOLD
120        I = KP2 - J
121        GDIF = OG(I) - OG(I-1)
122        TEMP2 = (G(I) - G(I-1)) - SIGMA*GDIF
123        TEMP3 = (C(I) - C(I-1)) + RMU*GDIF
124        DO 110 L = 1,NEQN
125          YOUT(L) = YOUT(L) + TEMP2*PHI(L,I)
126 110      YPOUT(L) = YPOUT(L) + TEMP3*PHI(L,I)
127 120    CONTINUE
128      DO 130 L = 1,NEQN
129        YOUT(L) = ((1.0 - SIGMA)*OY(L) + SIGMA*Y(L)) +
130     1             H*(YOUT(L) + (G(1) - SIGMA*OG(1))*PHI(L,1))
131 130    YPOUT(L) = HMU*(OY(L) - Y(L)) +
132     1                (YPOUT(L) + (C(1) + RMU*OG(1))*PHI(L,1))
133C
134      RETURN
135      END
136