function nonlinearBvpHW()

%   epsilon*u''(t) + u*(u' - 1) = 0     t in (0,1)
%
%         u(0) = alpha = -1
%         u(1) = beta  = 1.5
% 

         alpha = -1.0;
          beta =  1.5;
        epsilon = 0.02;  


   maxiter = 50;
         N = 200;
         x = linspace(0,1,N)';
         u = alpha + (beta - alpha)*x;     % initial quess   
         h = x(2) - x(1);
         h2 = h^2;
         f = zeros(N,1);
         
         
         tol = 10e-12;
         iter = 0;
         du = 1;
          I = 2:N-1;

   while ( norm(du,inf)> tol )
   
                f(1) = (epsilon/h2)*(alpha - 2*u(1) + u(2)) + (u(1).*u(2))/(2*h) - (u(1).*alpha)/(2*h) - alpha;
                f(I) = (epsilon/h2)*(u(I-1) - 2*u(I) + u(I+1)) + (u(I).*u(I+1))/(2*h) - (u(I).*u(I-1))/(2*h) - u(I);
                f(N) = (epsilon/h2)*(u(N-1) - 2*u(N) + beta) + (u(N).*beta)/(2*h) - (u(N).*u(N-1))/(2*h) - beta;
                
                 up = diag((epsilon/h2 + u(1:N-1)/(2*h)),1);
                 
                low = diag((epsilon/h2 - u(2:N)/(2*h)),-1);
                
                 d1 = -(2*epsilon)/h2 + u(2)/(2*h) - alpha/(2*h) - 1;
                 d2 = -(2*epsilon)/h2 + u(I+1)/(2*h) - u(I-1)/(2*h) - 1;
                 d3 = -(2*epsilon)/h2 + beta/(2*h) - u(N-1)/(2*h) - 1;     
                 
                  d = diag([d1 d2' d3]); 
                  
                  J = ( d + up + low );
                  
            du = -J\f;
             u = u + du;
             
          iter = iter+1;
          if iter==maxiter, disp('exceeded maximum iterations'), break;  end
          
          plot(x,u)
          pause
       
   end
   
   iter