function [A,b] = mkp(P,T,k,f,idxNeum,tauNeum,ID,uDirich)

nvlnka = size(P,2);
m = size(T,2);
mN = size(idxNeum,2);
Avlnka = zeros(nvlnka,nvlnka);
cvlnka = zeros(nvlnka,1);

for e=1:m
    idxe = T(:,e);
    xe = P(:,idxe);
    Ae = alok2d(xe,k(e));
    Avlnka(idxe,idxe) = Avlnka(idxe,idxe) + Ae;
    ce = clok2d(xe,f(e));
    cvlnka(idxe) = cvlnka(idxe) + ce;
end

% Neumannovy podminky
for s=1:mN
    idxs = idxNeum(:,s);
    xs = P(:,idxs);
    ds = dlok2d(xs,tauNeum(s));
    cvlnka(idxs) = cvlnka(idxs) + ds;
end

% Dirichletovy podminky
idxNonDirich = setdiff(1:nvlnka,ID);
A = Avlnka(idxNonDirich,idxNonDirich);
b = cvlnka(idxNonDirich);
b = b - Avlnka(idxNonDirich,ID)*uDirich;