function u = multigrid(l,r)

global Matice Interpolace Omegy

A = Matice{l};
n = length(r);
I = Interpolace{l};
omega = Omegy{l};

m = 3; % smoothing steps (of Richardon's method)
% [l,norm(r)]

if l==1
    
    % coarsest level correction
    u = A\r;
    
else

    u = zeros(n,1);
    
    % pre-smoothing
    for i=1:m
        u = u + omega*r;
        r = r - omega*A*r;
    end
%     [l,norm(r)]
    
    
    % coarse-grid correction
    rcoarse = I'*r;
    v = multigrid(l-1,rcoarse);
    u = u + I*v;
    r = r - A*I*v;
%     [l,norm(r)]
    
    % post-smoothing
    for i=1:m
        u = u + omega*r;
        r = r - omega*A*r;
    end
%     [l,norm(r)]
    
end