% twoPointBvpNeumannFD2
%
%   u_xx = e^x         on (0,1)
%
%      u'(0) = alpha
%       u(1) = beta
%
%  second order boundary treatment
%  2nd order one-sided difference    u_x = (-3*u_0 + 4*u_1 - u_2)/(h)

   clear, home, close all
    N = 40;
    x = linspace(0,1,N)';
   dm = dmSecond3Pt(x);
    h = x(2)-x(1);

   alpha = -1;       %  u'(0) = alpha
   beta = 0.2;       %  u (0) = beta
   
   uExact = exp(x) + (beta - exp(1) - alpha + 1) + x*(alpha - 1);
        f = exp(x);
   
   
    dm(1,1) = -1.5/h; dm(1,2) = 2/h; dm(1,3) = -0.5/h;   % Neumann BC at x=0
       f(1) = alpha;
       
    dm(N,:) = 0;                                 %  zero dirichlet bc at x=1
    dm(N,N) = 1;
       f(N) = beta;                     
     
      uh = dm\f;
      
      norm(uExact-uh,inf)
      plot(x,uh,'b--.',x,uExact,'r')