C FIGURE 5.8.1 
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                
      EXTERNAL FUNF,FUNA,FUNB
      CHARACTER*1 HORS                                                  
C                              N_ = NUMBER OF SUBINTERVALS           
C                              M_ = NUMBER OF UNKNOWNS
C                              NSTEP = NUMBER OF TIME STEPS
      PARAMETER (N_=10,M_=2*N_+1,NSTEP=200)
C                              BREAK POINTS ARE XPTS(K), K=0,N          
      COMMON /BLKX/ N,XPTS(0:N_)                                     
C                              UNKNOWNS ARE A(I), I=1,M                 
      COMMON /BLKA/ M,A(M_)
C                              PASS T(k+1) THROUGH COMMON
      COMMON /BLKT/ TKP1
      DIMENSION AMAT(M_,-3:3),F(M_)
      N = N_                                                         
      M = M_ 
      DO 5 K=0,N                                                        
         XPTS(K) = DBLE(K)/N
    5 CONTINUE                                                          
C                              INITIAL VALUES SET          
      DO 10 I=1,M                                                       
         CALL LOCATE(I,HORS,K)                                          
         IF (HORS.EQ.'H') A(I) =  EXP(-XPTS(K))
         IF (HORS.EQ.'S') A(I) = -EXP(-XPTS(K))
   10 CONTINUE                                                          
      DT = 1.D0/NSTEP
      DO 40 KT=1,NSTEP
         TKP1 = KT*DT
C                              BACKWARD DIFFERENCE METHOD REQUIRES
C                              SOLUTION OF (BMAT + DT*AMAT)*a(k+1) 
C                              = DT*b + BMAT*a(k) EACH TIME STEP
         L = 3
         DO 20 K=1,M                                                    
            F(K) = DT*GAUSS(FUNF,K,0)
            DO 15 I=MAX(K-L,1),MIN(K+L,M)
               AMATKI = GAUSS(FUNA,K,I) +
     &           EXP(0.5*TKP1)*PHI(K,0,1.D0)*PHI(I,0,1.D0)
               BMATKI = GAUSS(FUNB,K,I)
               F(K) = F(K) + BMATKI*A(I)
               AMAT(K,I-K) = BMATKI + DT*AMATKI
   15       CONTINUE                                                    
   20    CONTINUE
C                              SOLVE LINEAR SYSTEM USING BAND SOLVER    
         CALL LBAND(AMAT,A,F,M,L)
C                              CALCULATE MAXIMUM ERROR                  
         ERMAX = 0.0                                                    
         DO 30 J=0,100                                                  
            X = J/100.D0
            ERR = ABS(USOL(0,X) - EXP(TKP1-X))
            ERMAX = MAX(ERMAX,ERR)                                      
   30    CONTINUE                                                       
         PRINT 35, TKP1,ERMAX
   35    FORMAT (F8.3,E15.5)
   40 CONTINUE                                                          
      STOP                                                              
      END                                                               
      FUNCTION FUNF(X,K,I)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON /BLKT/ T
C                              INTEGRAND FOR K-TH ELEMENT OF "b" VECTOR
C                              (I IS NOT USED)
      OM  = OMEGA(0,X)
      OMt = OMEGA(0,X)
      OMx = OMEGA(1,X)
      PK  = PHI(K,0,X)
      PKx = PHI(K,1,X)
      FUNF = EXP(0.5*X**2*T)*(-OMt*PK - OMx*PKx + X*T*OM*PK)
      RETURN                                                            
      END                                                               
      FUNCTION FUNA(X,K,I)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON /BLKT/ T
C                              INTEGRAND FOR (K,I)-TH ELEMENT           
C                              OF "A" MATRIX                       
      PK  = PHI(K,0,X)
      PKx = PHI(K,1,X)
      PI  = PHI(I,0,X)
      PIx = PHI(I,1,X)
      FUNA = EXP(0.5*X**2*T)*(PKx*PIx - X*T*PK*PI)
      RETURN                                                            
      END                                                               
      FUNCTION FUNB(X,K,I)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON /BLKT/ T
C                              INTEGRAND FOR (K,I)-TH ELEMENT           
C                              OF "B" MATRIX                       
      PK  = PHI(K,0,X)
      PI  = PHI(I,0,X)
      FUNB = EXP(0.5*X**2*T)*PK*PI
      RETURN                                                            
      END                                                               
      SUBROUTINE LOCATE(I,HORS,K)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                
      CHARACTER*1 HORS                                                  
C                              DETERMINE THE TYPE ('H' OR 'S') OF BASIS 
C                              FUNCTION I, AND THE POINT (XPTS(K)) ON   
C                              WHICH IT IS CENTERED.                    
      IF (MOD(I,2).EQ.0) HORS = 'H'                                     
      IF (MOD(I,2).EQ.1) HORS = 'S'                                     
      K = I/2                                                           
      RETURN                                                            
      END                                                               
      FUNCTION OMEGA(IDER,X)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON /BLKT/ T
C                              EVALUATE OMEGA (IDER=0) OR ITS IDER-TH   
C                              DERIVATIVE AT X                          
      OMEGA = EXP(T)*HERM('H',0,IDER,X)
      RETURN                                                            
      END