function u = twogrid(r)

global A1 Lambda1 I1 A0 smoothingSteps

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 = A0\r0;
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
