C FIGURE 5.11.1
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                
      EXTERNAL FUNA,FUNB                                                
      CHARACTER*1 HORS                                                  
C                              N_ = NUMBER OF SUBINTERVALS           
C                              M_ = NUMBER OF UNKNOWNS
      PARAMETER (N_=4,M_=2*N_+1)                                       
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_)                                      
      DIMENSION AMAT(M_,-3:3),F(M_),AOLD(M_)
      EXACT(X) = EXP(X)*SIN(2.0287579*X)                                
      N = N_                                                         
      M = M_                                                     
      ALPHA = 4.0                                                       
      DO 5 K=0,N                                                        
         XPTS(K) = DBLE(K)/N                                            
    5 CONTINUE                                                          
C                              INITIAL VALUES ASSIGNED                  
      DO 10 I=1,M                                                       
         CALL LOCATE(I,HORS,K)                                          
         IF (HORS.EQ.'H') A(I) = 1.0                                    
         IF (HORS.EQ.'S') A(I) = 0.0                                    
   10 CONTINUE                                                          
C                              DO ITMAX INVERSE POWER ITERATIONS        
      ITMAX = 15                                                        
      DO 45 ITER=1,ITMAX                                                
C                              SAVE PREVIOUS A VECTOR                   
         DO 15 I=1,M                                                    
            AOLD(I) = A(I)                                              
   15    CONTINUE                                                       
C                              CALCULATE MATRIX A-ALPHA*B, AND F=B*AOLD 
         L = 3
         DO 25 K=1,M                                                    
            F(K) = 0.0                                                  
            DO 20 I=MAX(K-L,1),MIN(K+L,M)
               AMAT(K,I-K) = GAUSS(FUNA,K,I)
               BMAT = GAUSS(FUNB,K,I)
               F(K) = F(K) + BMAT*A(I)
               AMAT(K,I-K) = AMAT(K,I-K) - ALPHA*BMAT
   20       CONTINUE                                                    
   25    CONTINUE                                                       
C                              SOLVE LINEAR SYSTEM USING BAND SOLVER    
         CALL LBAND(AMAT,A,F,M,L)                                     
C                              ESTIMATE EIGENVALUE                      
         ANUM = 0.0                                                     
         ADEN = 0.0                                                     
         DO 30 I=1,M                                                    
            ANUM = ANUM + F(I)*AOLD(I)                                  
            ADEN = ADEN + F(I)*A(I)                                     
   30    CONTINUE                                                       
         EIGEN = ANUM/ADEN + ALPHA                                      
C                              CALCULATE MAXIMUM ERROR                  
         ERMAX = 0.0                                                    
         DO 35 J=0,100                                                  
            X = J/100.D0                                                
            ERR =
     &      ABS(USOL(0,X)/USOL(0,1.D0) - EXACT(X)/EXACT(1.D0))
            ERMAX = MAX(ERMAX,ERR)                                      
   35    CONTINUE                                                       
         PRINT 40, ITER,EIGEN,ERMAX                                     
   40    FORMAT (I5,2E17.8)                                             
   45 CONTINUE                                                          
      STOP                                                              
      END                                                               
      FUNCTION FUNA(X,K,I)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                
C                              INTEGRAND FOR (K,I)-TH ELEMENT OF        
C                              "A" MATRIX
      PK  = PHI(K,0,X)
      PKx = PHI(K,1,X)
      PI  = PHI(I,0,X)
      PIx = PHI(I,1,X)
      FUNA = EXP(-2.0*X)*(PKx*PIx-PK*PI)
      RETURN                                                            
      END                                                               
      FUNCTION FUNB(X,K,I)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                
C                              INTEGRAND FOR (K,I)-TH ELEMENT OF        
C                              "B" MATRIX
      PK  = PHI(K,0,X)
      PI  = PHI(I,0,X)
      FUNB = EXP(-2.0*X)*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)                                
C                              EVALUATE OMEGA (IDER=0) OR ITS IDER-TH   
C                              DERIVATIVE AT X                          
      OMEGA = 0.0                                                       
      RETURN                                                            
      END