function u = multigrid(l,r)

global A_array Lambda_array I_array A0 smoothingSteps

if l==0
    u = A0\r;
    return;
end

A1 = A_array{l};
Lambda1 = Lambda_array{l};
I1 = I_array{l};

n = length(r);
u = zeros(n,1);

% presmoothing
for i=1:smoothingSteps
    u = u + (1/Lambda1)*r;
    r = r - (1/Lambda1)*A1*r;
end

% coarse grid correction
r0 = I1'*r;
v0 = multigrid(l-1,r0);    % Zmena oproti twogrid
v1 = I1*v0;
u = u + v1;
r = r - A1*v1;

% postsmoothing
for i=1:smoothingSteps
    u = u + (1/Lambda1)*r;
    r = r - (1/Lambda1)*A1*r;
end
