C***********************************************************************
C
C   PROGRAM:  SAMPLE 1: STRIPPED DOWN PROGRAM USING CRANK-NICHOLSON
C             SUBROUTINES CNSET & CNSTEP() TO VALUE PUT OPTIONS ON
C             A NON-DIVIDEND PAYING SECURITY.
C
C***********************************************************************

      implicit double precision (a-h,k-l,o-z)
      parameter (N = 40)
      dimension U(0:N), ARR(0:N,4), PARM(15)

C**** SET MODEL PARAMETERS *********************************************

      open(4,file='samp1.out')

C     Crank-Nicholson algorithm parameters:

      data IMIN, IMAX, SMIN, SMAX, K / 0, 0, 0.0, 100.0, .025 /
      H = ( SMAX - SMIN ) /  N

C     Financial model parameters:

      data TMAT, R, SIGMA, XPRICE, IFUT /  .25, .10, .20, 50.0, 0 /
      PARM(1) = SIGMA
      PARM(2) = R

C**** SETUP TO SOLVE PDE ***********************************************

      call CNSET (N,SMIN,SMAX,K,IFN,IFUT,IMIN,IMAX,PARM,ARR)

      do I = 0, N
         S    = SMIN + I * H
         U(I) = max ( 0d0 , XPRICE - S )
      enddo
      
C**** SOLUTION LOOP TO SOLVE PDE ***************************************

      do 200 T = 0d0, TMAT-K/2, K
200      call CNSTEP ( T, U, ARR )

C**** PRINT RESULTS TO SCREEN AND FILE 4 *******************************

      do 300 I = 2, N-1, 2
         S    = SMIN + I * H
         Us   = (U(I+1) - U(I-1)) / (2 * H )
         Uss  = (U(I+1) - 2 * U(I) + U(I-1)) / ( H * H )
         write(4,350) S, U(I), Us, Uss
300      print 350, S, U(I), Us, Uss

350   format (f8.2, 3(10x, f10.4))

      stop
      end

C**** FUNCTION DEFINITIONS REQUIRED ************************************
C
C     The coefficients for the Black-Scholes option pricing model, with
C     S being the stock price, and R and D the assumed constant interest
C     rate and proportional dividend rate respectively, are 
C     FNA()  = SIGS * SIGS * S * S / 2.0
C     FNB()  = (R - D) * S
C     FNC()  = -R
C
C***********************************************************************

      double precision function COEFF()

      implicit double precision (a-h,k-l,o-z)
      dimension PARM(15)

      entry FNA(S,IFN,PARM)
         SIGMA  = PARM(1)
         FNA    = (SIGMA * S) ** 2 / 2.0
      return 
      
      entry FNB(S,IFN,PARM)
         R      = PARM(2)
         FNB    = R * S
      return

      entry FNC(S,IFN,PARM)
         FNC    = -R
      return
      
      entry FMIN(T,IFN,PARM)
         FMIN    = 0.0
      return
      
      entry FMAX(T,IFN,PARM)
         FMAX    = 0.0
      return
      end