function priklad2_1d_konvergenceMKP()

% Analyticke reseni:
% u(x) = 2*exp(-x^2)
%
% ulohy:
% -u''(x) = -2*exp(-x^2)*(4x^2-1) v (0,1)
% u(0) = 2
% 2*u'(1) = 2/e

uanal = @(x) (2*exp(-x.^2));
duanal = @(x) (-4*x.*exp(-x.^2));

x = linspace(0,1,100);
plot(x,uanal(x),'r');

for l=0:8
    m = 2*2^l;
    n = m+1; h = 1/m;
    poleh(l+1) = h;

    % diskretizace
    x = linspace(0,1,n);
    E = [1:(n-1);2:n];
    teziste = (x(1:(end-1))+x(2:end))/2;

    f = -2*exp(-teziste.^2).*(4*teziste.^2-2);
    k = ones(1,m);
    
    % okr. podm.
    Dirichlet = [1,n]; uD = [2,2/exp(1)];
    Neumann = []; tau = [];
    
    % MKP
    Ae = @(ke,xe) (ke/(xe(2)-xe(1))*[1,-1;-1,1]);
    be = @(fe,xe) (fe*(xe(2)-xe(1))/2*[1;1]);
    [A,b] = MKP(x,E,{Ae},{k},{be},{f},Dirichlet,uD,Neumann,tau);
    
    % reseni
    u = A\b;
    hold on
    plot(x,[uD(1),u',uD(2)])
    pause(2)

    u = [uD(1),u',uD(2)];
    teziste = (x(1:(end-1))+x(2:end))/2;
    L2chyba(l+1) = sqrt(h * sum((uanal(teziste)-0.5*u(1:(end-1))-0.5*u(2:end)).^2));

    H1chyba(l+1) = 0;
    for e=1:m
        xe = x([e,e+1]); ue = u([e,e+1]);
        duhe = (ue(2)-ue(1))/h;
        H1chyba(l+1) = H1chyba(l+1) + 0.5 * ( (duanal(xe(1))-duhe).^2 + (duanal(xe(2))-duhe).^2 );
    end
    H1chyba(l+1) = sqrt(h*H1chyba(l+1));
end
H1CelkovaChyba = sqrt(L2chyba.^2+H1chyba.^2);

figure
plot(log10(poleh),log10(L2chyba),'r',log10(poleh),log10(poleh.^2),'r--')
hold on
plot(log10(poleh),log10(H1chyba),'b',log10(poleh),log10(poleh.^1),'b--')
hold on
plot(log10(poleh),log10(H1CelkovaChyba),'g',log10(poleh),log10(poleh.^1),'g--')
