dx=.02; 
xArr=(0:dx:1); 
yArr=(0:dx:1); 

[xMesh,yMesh]=meshgrid(xArr,yArr); 
vMesh=0*xMesh; 
rhoMesh=0*xMesh; 

nIter=3000; 
for ii=1:nIter
    % apply boundary conditions
    vMesh(xMesh==0)=0; 
    vMesh(xMesh==1)=1; 
    
    vMesh(yMesh==0)=xMesh(yMesh==0);
    vMesh(yMesh==1)=xMesh(yMesh==0);
%     vMesh(yMesh==0)=0; 
%     vMesh(yMesh==1)=0; 

    vMesh(((xMesh-0.5).^2+(yMesh-0.5).^2)<=0.25^2)=0; 
    
    % iteration step
    for ix=2:length(xArr)-1
        for iy=2:length(yArr)-1
            vMesh(ix,iy)=(vMesh(ix+1,iy)+vMesh(ix-1,iy)+...
                vMesh(ix,iy+1)+vMesh(ix,iy-1)+rhoMesh(ix,iy)*dx^2)/4; 
        end 
    end 
    
    % plot
    if mod(ii,50)==0 || ii==nIter
        figure(1)
        clf
        surf(xMesh,yMesh,vMesh,'EdgeColor','none')
        set(gca,'FontSize',18)
        title(['iteration:',num2str(ii)])
        xlabel('x'); ylabel('y'); zlabel('V')
        colorbar
        view(2)
        drawnow 
    end 
end 

% save('tmp.mat','xMesh','yMesh','vMesh')