function [A,b] = MKP(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) + Ae(k(e),P(:,idx));
    end
    for i=1:nbe
        be = becka{i}; f = fka{i};
        bvlnka(idx) = bvlnka(idx) + be(f(e),P(:,idx));
    end
end
% Neumannovy okr. podm.
bvlnka(Neumann) = bvlnka(Neumann) + tau;
% 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
