function [A,b] = assemble(fce_tau,fce_q,fce_f,a,fa,x,e,I0)

% [A,b] = assemble(fce_tau,fce_q,fce_f,a,fa,x,e)
%          Sestavi matici (tuhosti) a vektor (zatizeni),
%          ktera vzejde z konecne-prvkove aproximace
%          nasledujici 1d variacni ulohy:
%           Najdi u z Vh : a(v,u) = f(v) pro vsechna v z Vh,
%          kde
%           a(v,u) = \int tau(x)*v'(x)*u'(x) dx +
%                    \int q(x)*v(x)*u(x) dx,
%           f(v) = \int f(x)*v(x) + fa*v(a),
%           Vh je dano Lagrangeovymi tvarovymi funkcemi
%           s nulovou okrajovou podminkou,
%           ktera je urcena indexy uzlu I0.

% pocet uzlu (pocet stupnu volnosti)
n = max(size(x));

% pocet elementu
ne = max(size(e));

A = zeros(n,n);
b = zeros(n,1);

for i=1:ne
    ei = e(i,:);
    x1 = x(ei(1));
    x2 = x(ei(2));
    
    taui = feval(fce_tau,x1+(x2-x1)/2);
    qi = feval(fce_q,x1+(x2-x1)/2);
    fi = feval(fce_f,x1+(x2-x1)/2);
    
    Aei = alok_tau(taui,x,ei) + alok_q(qi,x,ei);
    bei = blok_f(fi,x,ei) + blok_fa(a,fa,x,ei);
    
    nei = max(size(Aei));
    
    for j=1:nei
        for k=1:nei          
            if j==k | (~in(I0,ei(j)) & ~in(I0,ei(k)))
                A(ei(j),ei(k)) = A(ei(j),ei(k)) + Aei(j,k);
            end
        end
        if ~in(I0,ei(j))
            b(ei(j)) = b(ei(j)) + bei(j);
        end
    end
end
