SUBROUTINE INTRP(X, Y, XOUT, YOUT, YPOUT, NEQN, KOLD, PHI, PSI) * * The methods in the subroutine STEP approximate the solution near X by a * polynomial. Subroutine INTRP approximates the solution at XOUT by evaluating * the polynomial there. Information defining this polynomial is passed from * STEP so INTRP cannot be used alone. * * More to come!!! * IMPLICIT DOUBLE PRECISION (A-H, O-Z) DIMENSION Y(NEQN), YOUT(NEQN), YPOUT(NEQN), PHI(NEQN, 16) DIMENSION PSI(12) * DIMENSION G(13), W(13), RHO(13) DATA G(1) / 1.0D0 /, RHO(1) / 1.0D0 / * HI = XOUT - X KI = KOLD + 1 KIP1 = KI + 1 * * Initialize W(*) for computing G(*). * DO 5 I = 1, KI TEMP1 = DFLOAT(I) W(I) = 1.0D0 / TEMP1 5 CONTINUE TERM = 0.0D0 * * Compute G(*). * 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) TERM = PSIJM1 15 CONTINUE * * Interpolate. * DO 20 L = 1, NEQN YPOUT(L) = 0.0D0 YOUT(L) = 0.0D0 20 CONTINUE 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) YPOUT(L) = YPOUT(L) + TEMP3*PHI(L, I) 25 CONTINUE 30 CONTINUE DO 35 L = 1, NEQN 35 YOUT(L) = Y(L) + HI*YOUT(L) RETURN END