function [A,b] = FEM1d(l,k,f,tau)

n = length(f);
x = linspace(0,l,n+1);
Atilde = sparse(n+1,n+1);
btilde = zeros(n+1,1);

for e=1:n
    idxe = [e,e+1];
    he = x(idxe(2))-x(idxe(1));
    Ae = k(e)/he * [1,-1;-1,1];
    be = f(e)*he/2 * [1;1];
    Atilde(idxe,idxe) = Atilde(idxe,idxe) + Ae;
    btilde(idxe) = btilde(idxe) + be;
end
btilde(n+1) = btilde(n+1) + tau;
A = Atilde(2:n+1,2:n+1); b = btilde(2:n+1);
