function [ulevels,uexact] = asm1d(fce_tau,fce_q,fce_f,a,fa,xl,el,levels)

% function [ulevels,uexact] = asm1d(fce_tau,fce_q,fce_f,a,fa,xl,el,levels)

l = 0;
[Al,bl] = assemble(fce_tau,fce_q,fce_f,a,fa,xl,el,[1,max(size(xl))]);
ul = Al\bl;

for l=1:levels
    [xl,el,El] = refine(xl,el);

    ul = El*ul;

    [Al,bl] = assemble(fce_tau,fce_q,fce_f,a,fa,xl,el,[1,max(size(xl))]);
    
    Dl = diag(diag(Al)); % Jacobi smoother
    
    wl = inv(Dl)*bl;
    
    ul = ul + wl;
end

ulevels = ul;

% for test purposes only
uexact = Al\bl;
