function u = solve_poisson_multigrid(f)

% 1D Poisson equation with Dirichlet boundary conditions

N = length(f) - 1;      % N is a power of 2
h = 1/N;                % true grid: N+1 points
u = zeros(N+1,1);
I = 2:N;

if N > 2,
    
    % Run three iterations of Jacobi
    for k = 1:3
        u(I) = 1/2 * (u(I+1) + u(I-1) - h^2 * f(I));
    end
    
    % Recursive call at a lower scale to get an
    % approximation of the error
    residual = f(I) - 1/h^2 * (u(I+1) - 2*u(I) + u(I-1));
    residual = [0; residual; 0];
    coarse_residual = residual(1:2:end);
    coarse_err = solve_poisson_multigrid(-coarse_residual);
    err = zeros(N+1,1);
    err = interp1(1:2:(N+1), coarse_err, 1:(N+1), 'linear');
    u = u - err';
    
    % Three more iterations of Jacobi
    for k = 1:3
        u(I) = 1/2 * (u(I+1) + u(I-1) - h^2 * f(I));
    end    
       
else                    % N == 2
                        % 3 points
    
    u(I) = -1/2 * h^2 * f(I);
    
end

