function konvergenceMKP1d()

% -u''(x)) = 1 x z (0,1)
% u(0) = u(1) = 1
% analyt. reseni: 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)
    % MKP diskretizace
    h = hcka(i);
    x = 0:h:1;
    nvlnka = length(x);
    % data ulohy
    k = ones(1,nvlnka-1);
     f = ones(1,nvlnka-1);
    uDirich = [1;1]; idxDirich = [1,nvlnka];
    idxNonDirich = setdiff(1:nvlnka,idxDirich);
    tauNeum = []; idxNeum = [];
    % MKP
    [A,b] = mkp1d(x,k,f,idxNeum,tauNeum,idxDirich,uDirich);
    uvlnka = zeros(nvlnka,1);
    u = A\b;
    uvlnka(idxNonDirich) = u;
    uvlnka(idxDirich) = uDirich;
    du = (uvlnka(2:end)-uvlnka(1:end-1))'/h;
    xSimpson = [x(1:end-1);0.5*(x(1:end-1)+x(2:end));x(2:end)];
    wSimpson = [1/6,2/3,1/6];
    duanalyt = 0.5-xSimpson;
   err(i) = sqrt(sum(h*wSimpson*([du;du;du]-duanalyt).^2));
end

plot(log(hcka),log(err),'b',log(hcka),log(hcka),'r--')
