N = 100
nx = N-1;
U = zeros(nx+1); 
x = linspace(0,1,nx+1);
y = x;

% Boundary conditions
U(1,round(nx/4):round(nx/2))=1;
U(round(nx/2):round(3*nx/4),1)=2;

% Iterations
errest = 1
k = 1
while(errest>1.e-4)
    errest = 0;
    k = k+1
    for j=2:nx
        for i=2:nx
            Unew = 0.25*(U(i+1,j)+U(i-1,j)+U(i,j+1)+U(i,j-1));
            errest = max(abs(Unew-U(i,j)),errest);
            U(i,j) = Unew;
        end
    end
    if k>10000 
        break;
    end
    errest
end

% Plot the solution
surf(x,y,U)
shading('interp')
view(0,90)
axis('square')
colorbar()
    


