% bvp1dAdvDifFD2
% epsilon*u_xx - u_x = -1

   clear, home, close all
        N = 25;
   epsilon = 0.01;
     alpha = 1;      %  u(0) = alpha
      beta = 3;       %  u(1) = beta
   
     x = linspace(0,1,N)';
   dm2 = dmSecond3Pt(x);
   dm1 = dmFirst3Pt(x) ;
    dm = epsilon.*dm2 - dm1;

      u = alpha + x + (beta-alpha-1).*(exp(x./epsilon)-1)./(exp(1./epsilon)-1);   % exact solution
      
        f = -1*ones(N,1); 
      f(1) = alpha;
      f(N) = beta;
      dm(1,:) = 0;  dm(1,1) = 1;
      dm(N,:) = 0;  dm(N,N) = 1;
      
      uh = dm\f;
   
      norm(u-uh,inf)
      plot(x,uh,'b--o',x,u,'r')