function fin_diff_solve_neumann

% compute finite difference solution to u''+u=0 on uniform grid
% with u'(-1)=1 and u(1)=1
% exact solution is u(x)=(cos(x)+sin(x))/(cos(1)+sin(1))

% we repeat with N=4,8,16,...,256 to show convergence using two different
% methods for applying the Neumann boundary condition

m=7;
err1=zeros(1,m); err2=zeros(1,m);
for k=1:m
    N=2^(k+1);
    % set up uniform grid
    h=2/N;
    x=-1:h:1; x=x(:);

    % set up finite difference matrix
    e=ones(N+1,1);
    A=spdiags([e -2*e e], -1:1, N+1, N+1);
    A=A/h^2+speye(N+1);
    
    % apply Dirichlet boundary condition at x=1
    A(N+1,:)=[zeros(1,N),1];

    % set up right-hand-side (including Dirichlet boundary data)
    b=zeros(N-1,1);  b(N+1)=1;
    
    % apply Neumann boundary condition in form (U_1-U_0)/h=tan(1)
    A(1,:)=[-1/h,1/h,zeros(1,N-1)];
    b(1)=1;

    % solve linear system
    U=A\b;

    % compute error
    err1(k)=norm(U-(cos(x)+sin(x))/(cos(1)+sin(1)),Inf);

    % now apply Neumann boundary condition using fictitious node
    A(1,:)=[1-2/h^2,2/h^2,zeros(1,N-1)];
    b(1)=2/h;
    
    % solve linear system
    U=A\b;

    % compute error
    err2(k)=norm(U-(cos(x)+sin(x))/(cos(1)+sin(1)),Inf);
end

Ns=2.^(2:m+1);
loglog(Ns,err1,'bx-',Ns,err2,'rx-')
hold on
loglog(Ns,1./Ns,'b--',Ns,1./Ns.^2,'r--')
xlabel('N')
ylabel('Error')
legend('Convergence of FD scheme with 1st order BC','Convergence of FD scheme with 2nd order BC','O(h)','O(h^2)')