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

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

rhoFMesh=0*xMesh; 
epsMesh=1+0*xMesh; 

epsMesh(((xMesh-0.5).^2+(yMesh-0.5).^2)<=0.25^2)=16; 
rhoFMesh(((xMesh-0.5).^2+(yMesh-0.5).^2)<=0.05^2)=1; 

nIter=3000; 
for ii=1:nIter
    % apply boundary conditions
    vMesh(xMesh==0)=0; 
    vMesh(xMesh==1)=0; 
    
%     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)=(...
                (epsMesh(ix+1,iy)+epsMesh(ix,iy))*vMesh(ix+1,iy)+...
                (epsMesh(ix-1,iy)+epsMesh(ix,iy))*vMesh(ix-1,iy)+...
                (epsMesh(ix,iy+1)+epsMesh(ix,iy))*vMesh(ix,iy+1)+...
                (epsMesh(ix,iy)+epsMesh(ix,iy-1))*vMesh(ix,iy-1)+...
                rhoFMesh(ix,iy)*dx^2)/...
                (epsMesh(ix+1,iy)+epsMesh(ix-1,iy)+epsMesh(ix,iy+1)+...
                epsMesh(ix,iy-1)+4*epsMesh(ix,iy)); 
            exMesh(ix,iy)=-(vMesh(ix,iy+1)-vMesh(ix,iy-1))/2/dx; 
            eyMesh(ix,iy)=-(vMesh(ix+1,iy)-vMesh(ix-1,iy))/2/dx; 
        end 
    end 
    
    % plot
    if mod(ii,100)==0 || ii==nIter
        
        figure(1)
        clf
        contourf(xMesh,yMesh,vMesh,50,'EdgeColor','none')
        hold on 
        quiver(xMesh,yMesh,exMesh,eyMesh,2, 'color','w')
        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')