C FIGURE 5.3.3a
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                
      EXTERNAL FUNJ,FUNF                                                
      CHARACTER*1 HORS                                                  
C                              N_ = NUMBER OF SUBINTERVALS           
C                              M_ = NUMBER OF UNKNOWNS
      PARAMETER (N_=4,M_=2*N_)                                          
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 AJACOB(M_,-3:3),F(M_),DELTA(M_)
      DATA PI/3.14159265358979312D0/                                    
      N = N_                                                         
      M = M_                                                       
      DO 5 K=0,N                                                        
         XPTS(K) = K*PI/N                                               
    5 CONTINUE                                                          
C                              INITIAL GUESS FOR NEWTON METHOD          
      DO 10 I=1,M                                                       
         CALL LOCATE(I,HORS,K)                                          
         IF (HORS.EQ.'H') A(I) = 0.5 - XPTS(K)/PI                       
         IF (HORS.EQ.'S') A(I) = -1.0/PI                                
   10 CONTINUE                                                          
C                              DO ITMAX NEWTON ITERATIONS               
      ITMAX = 15                                                        
      DO 40 ITER=1,ITMAX                                                
C                              CALCULATE F VECTOR AND JACOBIAN          
         L = 3
         DO 20 K=1,M                                                    
            F(K) = GAUSS(FUNF,K,0)                                      
            DO 15 I=MAX(K-L,1),MIN(K+L,M)
               AJACOB(K,I-K) = GAUSS(FUNJ,K,I)
   15       CONTINUE                                                    
   20    CONTINUE                                                       
C                              SOLVE LINEAR SYSTEM USING BAND SOLVER    
         CALL LBAND(AJACOB,DELTA,F,M,L)                                   
C                              UPDATE SOLUTION VECTOR                   
         DO 25 I=1,M                                                    
            A(I) = A(I)-DELTA(I)                                            
   25    CONTINUE                                                       
C                              CALCULATE MAXIMUM ERROR                  
         ERMAX = 0.0                                                    
         DO 30 J=0,100                                                  
            X = J*PI/100.D0                                             
            ERR = ABS(USOL(0,X) - SIN(X+PI/6.D0))
            ERMAX = MAX(ERMAX,ERR)                                      
   30    CONTINUE                                                       
         PRINT 35, ITER,ERMAX                                           
   35    FORMAT (I5,E15.5)                                              
   40 CONTINUE                                                          
      STOP                                                              
      END                                                               
      FUNCTION FUNF(X,K,I)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                
C                              INTEGRAND FOR K-TH ELEMENT OF F VECTOR
C                              (I IS NOT USED)
      U  = USOL(0,X)
      Ux = USOL(1,X)
      PK  = PHI(K,0,X)
      PKx = PHI(K,1,X)
      FUNF = -Ux*PKx -Ux**2*PK - U**2*PK + U*PK + PK
      RETURN                                                            
      END                                                               
      FUNCTION FUNJ(X,K,I)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                
C                              INTEGRAND FOR (K,I)-TH ELEMENT           
C                              OF JACOBIAN MATRIX                       
      U  = USOL(0,X)
      Ux = USOL(1,X)
      PK  = PHI(K,0,X)
      PKx = PHI(K,1,X)
      PI  = PHI(I,0,X)
      PIx = PHI(I,1,X)
      FUNJ = -PKx*PIx - 2*Ux*PK*PIx - 2*U*PK*PI + 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 = 'S'                                     
      IF (MOD(I,2).EQ.1) HORS = 'H'                                     
      IF (I.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 = 0.5*HERM('H',0,IDER,X) - 0.5*HERM('H',N,IDER,X)           
      RETURN                                                            
      END