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

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

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
    
    for s=1:smoothingsteps
        dls = bl - Al*ul;
        wls = inv(Dl)*dls;
        ul = ul + wls;
    end 
end

ulevels = ul;

% for test purposes only
uexact = Al\bl;
