create multigrid solver for 2d poissons equation with gauss seidel method used. in matlab

main.m
function u = multigrid_solver_poisson(f, u0, N, max_levels, nu1, nu2)
    h = 1/(N+1);
    x = h:h:1-h;
    
    u = cell(max_levels,1);
    for l = max_levels:-1:2
        Nl = N/(2^(l-1));
        hl = 1/(Nl+1);
        
        if l == max_levels
            u{l} = GS(u0, f, N, nu1);
        else
            u{l} = GS(u{l+1}, f, N, nu1);
        end
        
        for i = 1:nu2
            u{l} = GS(u{l}, f, N, 1);
        end
        
        u{l} = prolong(u{l+1});
        difference = u{l} - u{l};
        
        for i = 1:nu2
            u{l} = GS(u{l}, f, N, 1);
        end
        
        u{l} = u{l} + difference;
    end
    
    u{1} = GS(u{2}, f, N, nu1);
end

function u = GS(u, f, N, nu)
    h = 1/(N+1);
    x = h:h:1-h;
    
    for iter = 1:nu
        for i = 2:N+1
            for j = 2:N+1
                u(i,j) = (f(i,j)*h^2 + u(i-1,j) + u(i+1,j) + u(i,j-1) + u(i,j+1))/4;
            end
        end
    end
end

function u = prolong(u)
    N = size(u,1) - 2;
    Nl = size(u,1) - 2;

    unew = zeros(2*N+1, 2*N+1);
    unew(1:2:end, 1:2:end) = u;
    unew(2:2:end-1, 1:2:end) = 0.5*(u(1:end-1,:) + u(2:end,:));
    unew(1:2:end, 2:2:end-1) = 0.5*(u(:,1:end-1) + u(:,2:end));
    unew(2:2:end-1, 2:2:end-1) = 0.25*(u(1:end-1,1:end-1) + u(2:end,1:end-1) + u(1:end-1,2:end) + u(2:end,2:end));
    u = unew;
end
1344 chars
57 lines