SUBROUTINE CNSET (IN,SMIN,SMAX,K,IFN,IFUT,ISMN,ISMX,PARM,ARR)

C***********************************************************************
C
C  SUBROUTINE CNSET (...)
C
C  Subroutine sets up coefficient array needed for Crank-Nicholson
C  algorithm to solve 1 state variable plus time partial differential
C  equations.  Use in conjunction with routine CNSTEP.  PDE has form  
C
C                FNA * Uss + FNB * Us + FNC * U - Ut = 0
C
C  Arguments: IN      number of grid intervals in state space S 
C             SMIN    minimum value of state variable S
C             SMAX    maximum value of state variable S
C             K       step size in time direction
C             IFN     flag available for passing to coeff. fcns.
C             IFUT    flag setting FNC = 0 for futures contract pricing
C             ISMN    flag for SMIN boundary (0 for quadratic extrapol.)
C             ISMX    flag for SMAX boundary (1 for given values       )
C             PARM    vector of model parameters for coeff. fcns.
C             ARR     output array of coefficients for CNSTEP
C                     dimension (4,IN+1) by calling program
C
C  Other routines called:  functions FNA, FNB, FNC(S,IFN,PARM) must be
C                     externally defined and available to subroutine.
C
C  Author:     R. A. Jones     15 December 1988
C
C***********************************************************************
 
      IMPLICIT  DOUBLE PRECISION ( A-H, K-L, O-Z )
      COMMON   /CNCOM/ N,ISMIN,ISMAX,IIFN,XPARM(15)
      DIMENSION PARM( 15 ), ARR( 0:IN, 4 )

      N      = IN 
      IIFN   = IFN
      ISMIN  = ISMN
      ISMAX  = ISMX
      DO 50 I = 1, 15
  50     XPARM(I) = PARM(I)
      H      = ( SMAX - SMIN ) / DBLE( N )
      FUTURE = 1D0
      IF ( IFUT .EQ. 1)    FUTURE = 0D0

C**** FIRST DO 'INTERIOR' COEFFICIENTS ********************************

      DO 100  I = 1, N-1

         S        =  SMIN +  DBLE(I) * H
         AX       =  FNA(S,IFN,PARM) * 2D0 * K
         BX       =  FNB(S,IFN,PARM) * H * K
         CX       =  FNC(S,IFN,PARM) * FUTURE * 2D0 * H * H * K
         DENOM    =  CX - 2D0 * AX - 4D0 * H * H

         ARR(I,1) =  ( AX - BX ) / DENOM
         ARR(I,2) =  1D0
         ARR(I,3) =  ( AX + BX ) / DENOM
         ARR(I,4) =  1D0 + 8D0 * H * H / DENOM
                      
 100  CONTINUE

C**** THEN HANDLE BOUNDARIES ACCORDING TO FLAGS ************************

      IF ( ISMIN .EQ. 1 )  THEN

C        CASE OF KNOWN VALUE AT SMIN: ISMIN = 1
         ARR(0,1) =  0D0
         ARR(0,2) =  1D0
         ARR(0,3) =  0D0

      ELSE

C        CASE OF QUADRATIC EXTRAPOLATION AT SMIN: ISMIN = 0
         G        =  ARR(1,3) / ( ARR(2,2) + 3D0 * ARR(2,3) )
         ARR(0,1) =  0D0
         ARR(0,2) =  G * ARR(2,3) - ARR(1,1)
         ARR(0,3) =  G * ( ARR(2,1) - 3D0 * ARR(2,3) ) - ARR(1,2)
         ARR(0,4) =  G

      ENDIF

      IF ( ISMAX .EQ. 1 )  THEN

C        CASE OF KNOWN VALUE AT SMAX: ISMAX = 1
         ARR(N,1) =  0D0
         ARR(N,2) =  1D0
         ARR(N,3) =  0D0

      ELSE

C        CASE OF QUADRATIC EXTRAPOLATION AT SMAX: ISMAX = 0
         G        =  ARR(N-1,1) / ( ARR(N-2,2) + 3D0 * ARR(N-2,1) )
         ARR(N,1) =  G * ( ARR(N-2,3) - 3D0 * ARR(N-2,1) ) - ARR(N-1,2)
         ARR(N,2) =  G * ARR(N-2,1) - ARR(N-1,3)
         ARR(N,3) =  0D0
         ARR(N,4) =  G

      ENDIF

      RETURN
      END
      
      SUBROUTINE CNSTEP ( T, U, ARR )

C***********************************************************************
C
C  SUBROUTINE CNSTEP (...)
C  
C  Subroutine takes 1 step in time direction in solving 1 state variable
C  PDE using Crank-Nicholson algorithm.  T is current time used only for
C  passing to boundary value functions FMIN(T) and FMAX(T) if ISMIN or
C  ISMAX are set to 1.  U(0:N) is N+1 dimensional vector of solution so
C  far.  ARR() is coefficient array set up by CNSET().
C
C***********************************************************************

      IMPLICIT   DOUBLE PRECISION ( A-H, K-L, O-Z )
      COMMON    /CNCOM/ N, ISMIN, ISMAX, IFN, PARM(15)
      DIMENSION  ARR( 0:N, 4 ), U( 0:N )

C     NOTE: PARAMETER NMAX MUST BE .GE. N FOR TRIDAG ALGORITHM
      PARAMETER ( NMAX = 1000 )
      COMMON    /TRICOM/ D( 0:NMAX ), GAM( 0:NMAX )

C     SET UP RIGHT HAND SIDE OF SYSTEM TRIDIAGONAL SYSTEM (ABC)U = D

      DO 100  I = 1, N-1
         D(I) = - ARR(I,1)*U(I-1) - ARR(I,4)*U(I) - ARR(I,3)*U(I+1)
 100  CONTINUE

      IF ( ISMIN .EQ. 1 )  THEN
C        GET SOLUTION VALUE AT RMIN
         D(0) = FMIN(T,IFN,PARM)
      ELSE
         D(0) = D(2) * ARR(0,4) - D(1)
      ENDIF

      IF ( ISMAX .EQ. 1 )  THEN
C        GET SOLUTION VALUE AT RMAX
         D(N) = FMAX(T,IFN,PARM)
      ELSE
         D(N) = D(N-2) * ARR(N,4) - D(N-1)
      ENDIF

      CALL TRIDAG ( ARR(0,1), ARR(0,2), ARR(0,3), D, GAM, U, N )

      RETURN
      END
      
C**** TRIDIAGONAL SOLN. ALGORITHM FROM "NUMERICAL RECIPES", P. 40 ******
 
C     SOLVES: (ABC)X = D FOR X.  N=DIMENSION.  A,B,C,D, NOT ALTERED
C     NOTE: SUBSCRIPTS RUN FROM 0 AND SCRATCH VECTOR GAM VARIABLE DIMEN.
 
      SUBROUTINE TRIDAG ( A, B, C, D, GAM, X, N )

      IMPLICIT  DOUBLE PRECISION ( A-H, K-L, O-Z )
      DIMENSION A(0:*), B(0:*), C(0:*), D(0:*), GAM(0:*), X(0:*)

      BET  = B(0)
      X(0) = D(0) / BET
      DO 10 J = 1, N
         GAM(J) = C(J-1) / BET
         BET    =  B(J) - A(J) * GAM(J)
         X(J)   = (D(J) - A(J) * X(J-1)) / BET
  10  CONTINUE

      DO 20 J = N-1, 0, -1
  20     X(J)   =  X(J) - GAM(J+1) * X(J+1)
      RETURN
      END