function convergenceFEM1d()

% -u''(x)) = 1 x z (0,1)
% u(0) = u(1) = 1
% analyt. solution: u(x) = 1-0.5*x*(x-1)

xx = 0:0.001:1;
uanalyt = 1-0.5*xx.*(xx-1);
plot(xx,uanalyt,'r');

hcka = 0.5*2.^[0:-1:-5];

for i=1:length(hcka)
    % FEM solution
    h = hcka(i);
    P = 0:h:1;
    n = length(P); m = n-1;
    E = [1:m;2:m+1];
    ID = [1,m+1]; uhat = [1,1];
    IN = []; tauhat = [];
    k = ones(1,m);
    f = ones(1,m);
    [A,b,utilde] = FEM(P,E,{'aloc_1dLaplace'},{k},{'bloc_1dsource'},{f},ID,uhat,IN,tauhat);
    uH = A\b;
    u = utilde;
    nonDirNodes = setdiff(1:m+1,ID);
    u(nonDirNodes) = u(nonDirNodes) + uH;
    hold on
    plot(P,u')
    % FEM errors in H1-seminorm and L2-norm
    du = (u(2:end)-u(1:end-1))'/h;
    xTrapezoidal = [P(1:end-1);P(2:end)];
    wTrapezoidal = [1/2,1/2];
    duanalyt = 0.5-xTrapezoidal;
    H1err(i) = sqrt(sum(h*wTrapezoidal*([du;du]-duanalyt).^2));
    xSimpson = [P(1:end-1);0.5*(P(1:end-1)+P(2:end));P(2:end)];
    wSimpson = [1/6,2/3,1/6];
    uanalyt = 1-0.5*xSimpson.*(xSimpson-1);
    L2err(i) = sqrt(sum(h*wSimpson*([u(1:end-1)';0.5*(u(1:end-1)+u(2:end))';u(2:end)']-uanalyt).^2));

    disp(sprintf('h = %1.5g: H1err = %1.5g, L2err = %1.5g',h,H1err(i),L2err(i)))
    pause
end

figure
plot(log(hcka),log(H1err),'b',log(hcka),log(hcka),'b--',log(hcka),log(L2err),'r',log(hcka),log(hcka.^2),'r--')
