function [A,b] = MKP_2d(P,E,Aecka,kcka,becka,fka,Dirichlet,uD,Neumann,tau)

d = size(P,1);
nAe = length(Aecka); nbe = length(becka);
n = size(P,2); m = size(E,2);
Avlnka = zeros(n,n);
bvlnka = zeros(n,1);
for e=1:m
    idx = E(:,e);
    for i=1:nAe
        Ae = Aecka{i}; k = kcka{i};
        Avlnka(idx,idx) = Avlnka(idx,idx) + feval(Ae,k(e),P(:,idx));
    end
    for i=1:nbe
        be = becka{i}; f = fka{i};
        bvlnka(idx) = bvlnka(idx) + feval(be,f(e),P(:,idx));
    end
end

% Neumannovy okr. podm.
mN = size(Neumann,2);
for e=1:mN
    INe = Neumann(:,e);
    Pe = P(:,INe);
    de = tau(e)*norm(Pe(:,2)-Pe(:,1))/2;
    bvlnka(INe) = bvlnka(INe) + de;
end

% Dirichletovy okr. podm.
neDirichlet = setdiff(1:n,Dirichlet);
A = Avlnka(neDirichlet,neDirichlet); b = bvlnka(neDirichlet);
if length(Dirichlet)>0
    Astriska = Avlnka(neDirichlet,Dirichlet); b = b - Astriska*uD';
end
