function u = solvepoissonmultigrid(f)

% 1D Poisson equation with Dirichlet boundary conditions

N = length(f) - 1;     % N should be a power of 2
h = 1/N;                  % complete grid with endpoints: N+1 points
u = zeros(N+1,1);     % Initial guess
I = 2:N;                    % Interior indices
omega = 1/4;           % Relaxation parameter

if N > 2,
    
    % Weighted Jacobi iterations
    for k = 1:5
        u(I) = u(I) - omega * (-u(I+1)+2*u(I)-u(I-1) + h^2*f(I));
    end
    
    res = [0; f(I) - 1/h^2 * (u(I+1) - 2*u(I) + u(I-1)); 0];
    %res = [0; 1/4*res(I+1) + 1/2*res(I) + 1/4*res(I-1); 0];
    coarse_res = res(1:2:end);
    coarse_err = solve_poisson_multigrid(coarse_res);
    err = interp1(1:2:(N+1), coarse_err, 1:(N+1), 'linear');
    %err = [0, 1/4*err(I+1) + 1/2*err(I) + 1/4*err(I-1), 0];
    u = u + err';

    % Weighted Jacobi iterations
    for k = 1:5
        u(I) = u(I) - omega * (-u(I+1)+2*u(I)-u(I-1) + h^2*f(I));
    end
       
else                    % N == 2
                          % 3 points
    
    u(I) = -1/2 * h^2 * f(I);
    
end
