C FIGURE 5.4.2
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                
C                              NSUBS = NUMBER OF SUBINTERVALS           
      PARAMETER (NSUBS=4)                                               
C                              BREAK POINTS ARE XPTS(K), K=0,N          
      COMMON /BLKX/ N,XPTS(0:NSUBS)                                     
C                              UNKNOWNS ARE A(I), I=1,M                 
      COMMON /BLKA/ M,A(2*NSUBS)                                        
      DIMENSION AMAT(2*NSUBS,-2:2),B(2*NSUBS),BETA(2)
      N = NSUBS                                                         
      M = 2*NSUBS                                                       
      ALPHA = 3.636                                                     
      DO 5 K=0,N                                                        
         XPTS(K) = (DBLE(K)/N)**ALPHA                                   
    5 CONTINUE                                                          
C                              CALCULATE RIGHT HAND SIDE VECTOR AND     
C                              COEFFICIENT MATRIX                       
      BETA(1) = 0.5 - 0.5/SQRT(3.D0)                                    
      BETA(2) = 0.5 + 0.5/SQRT(3.D0)                                    
      L = 2 
      DO 15 LL=0,N-1                                                     
      DO 15 J=1,2                                                       
         K = 2*LL+J                                                      
         ZK = XPTS(LL) + BETA(J)*(XPTS(LL+1)-XPTS(LL))                     
         B(K) = -OMEGA(2,ZK) + 0.11D0*OMEGA(0,ZK)/ZK**2                 
         DO 10 I=MAX(K-L,1),MIN(K+L,M)
            AMAT(K,I-K) = PHI(I,2,ZK) - 0.11D0*PHI(I,0,ZK)/ZK**2
   10    CONTINUE                                                       
   15 CONTINUE                                                          
C                              SOLVE LINEAR SYSTEM USING BAND SOLVER    
      CALL LBAND(AMAT,A,B,M,L)                                        
C                              CALCULATE MAXIMUM ERROR                  
      ERMAX = 0.0                                                       
      DO 20 J=0,100                                                     
         X = J/100.D0                                                   
         ERR = ABS(USOL(0,X,1,M) - X**1.1D0)
         ERMAX = MAX(ERMAX,ERR)                                         
   20 CONTINUE                                                          
      PRINT 25, ERMAX                                                   
   25 FORMAT (E15.5)                                                    
      STOP                                                              
      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 /BLKX/ N,XPTS(0:1)                                         
C                              EVALUATE OMEGA (IDER=0) OR ITS IDER-TH   
C                              DERIVATIVE AT X                          
      OMEGA = 1.1D0*HERM('S',N,IDER,X)                                  
      RETURN                                                            
      END                                                               

      include "fig533b.f"
      include "fig044.f"