C FIGURE 4.2.1
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                
C                             N = NUMBER OF GRID POINTS                 
      PARAMETER (N=100)                                                 
      DIMENSION A(N),B(N),C(N),U(0:N),DELTA(N),F(N)                     
      DATA PI/3.14159265358979312D0/                                    
      DX = PI/N                                                         
C                             BOUNDARY VALUES SET                       
      U(0) =  0.5                                                       
      U(N) = -0.5                                                       
C                             INITIAL GUESS SATISFIES BNDRY. CONDITIONS 
      DO 10 I=1,N-1                                                     
         XI = I*DX                                                      
         U(I) = 0.5-XI/PI                                               
   10 CONTINUE                                                          
C                             BEGIN NEWTON ITERATION                    
      DO 50 ITER=1,15                                                   
         DO 20 I=1,N-1                                                  
C                             CALCULATE F AND 3 DIAGONALS OF JACOBIAN   
            F(I) = (U(I+1)-2*U(I)+U(I-1))/DX**2 -                       
     &            ((U(I+1)-U(I-1))/(2*DX))**2 - U(I)**2 + U(I) + 1      
            A(I) =  1/DX**2 + (U(I+1)-U(I-1))/(2*DX**2)                 
            B(I) = -2/DX**2 - 2*U(I) + 1                                
            C(I) =  1/DX**2 - (U(I+1)-U(I-1))/(2*DX**2)                 
   20    CONTINUE                                                       
C                             CALL TRI TO SOLVE TRIDIAGONAL SYSTEM      
         CALL TRI(A,B,C,DELTA,F,N-1)                                    
C                             UPDATE SOLUTION AND CALCULATE MAX. ERROR  
         ERMAX = 0.0                                                    
         DO 30 I=1,N-1                                                  
            U(I) = U(I) - DELTA(I)                                      
            XI = I*DX                                                   
            ERR = ABS(U(I) - SIN(XI+PI/6.))                             
            ERMAX = MAX(ERMAX,ERR)                                      
   30    CONTINUE                                                       
         PRINT 40, ITER,ERMAX                                           
   40    FORMAT (I5,E15.5)                                              
   50 CONTINUE                                                          
      STOP                                                              
      END