function [u, xx, yy] = myfd2dpoisson(rhs, bound, bsize, np, nit)
    h =  bsize/(np+1);
    h2 = h^2;
    xx = 0:h:bsize;
    yy = xx;

    u = zeros(np+2, np+2);  
    
    % Boundary condition
    for i = 1:np+2
        x = xx(i);
        y = yy(i);
        u(i, 1) = bound(x, yy(1));           % bottom
        u(i, np+2) = bound(x, yy(np+2));     % top
        u(1, i) = bound(xx(1), y);           % left
        u(np+2, i) = bound(xx(np+2), y);     % right
    end

    idxy = 2:(np+1);
    
    b = zeros(np+2, np+2);
    for i = idxy
        x = xx(i);
        for j = idxy
            y = yy(j);
            b(i, j) = h2*rhs(x, y);  %internal point only
        end
    end
    
    % Iteration loop
    for it = 1:nit
       u(idxy, idxy) = (u(idxy+1, idxy) + u(idxy-1, idxy) + u(idxy, idxy+1) + u(idxy, idxy-1) - b(idxy, idxy))/4;
    end
end
