% Boundary value problem
% Steady-state heat equation in 1 D
% Neumann boundary conditions
%
clear all
nsteps=4;
hh=zeros(nsteps,1);
a=rand;
for l=1:nsteps
    % for an increasing number of grid points m
    m=2^(3+l);
    e=ones(m+2,1);
    h=1/(m+1);
    hh(l)=h;
    % Assembling the matrix A as a sparse matrix
    A=spdiags([e -2*e e], -1:1,m+2,m+2);
    A(1,1)=-h;
    A(1,2)=h;
    %A(1,3)=h/2;
    A(end,end-1)=0;
    A(end,end)=h^2;
    A=1/(h^2)*A;
    % BCs
    alpha=sin(0)*pi*a;
    beta=-cos(pi*a);
    % rhs function
    f=inline('cos(pi.*g.*x).*pi^2.*g^2','x','g');
    % grid including boundaries
    xx=[0:h:1];
    % internal nodes, peeling off the boundaries
    xin=xx(2:end-1);
    %
    % rhs of the system conicides with the values of the rhs function
    % on the internal nodes of the grid
    %
    F=f(xx',a);
    %
    % Neumann BCs on the left boundary point implies:
    %
    F(1,1)=alpha+h/2*f(0,a);
    % Dirichlet BC on the right boundary point implies:
    F(end,1)=beta;
    % Solving the linear system
    U=A\F;
    plot(xx,U)
    hold on
    Uex=-cos(pi*a*xx');
    plot(xx,Uex,'dr')
    hold off
    %
    % Computing the 2-norm of the error as a piecewise constant function
    %
    err(l)=norm(U-Uex)*h^(1/2);
    %
end
loglog(hh,err,'s-')
hold on
loglog(hh,hh.^2,'k')
%err'./(hh.^2)
%
xlabel('step-size')
ylabel('2-norm of the error')
title('Log-log plot of the 2-norm of the error for a test BVP problem, Neumann BCs: order 2 discretization')
hold off