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

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

nIter=4000; 
beta=40; 

for ii=1:nIter
    % apply boundary conditions
    vMesh(yMesh==0)=0; 
    vMesh(yMesh>=xMesh*tand(beta))=0; 
    vMesh(xMesh==1)=1;

    % 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 

%% post-plot
rArr=10.^(-2:0.05:0); phi0=beta/4; 
vArr=interp2(xMesh,yMesh,vMesh,rArr*cosd(phi0),rArr*sind(phi0));
vAnal=rArr.^(180/beta); 

figure(30)
loglog(rArr,vArr, '-o', rArr,vAnal)
set(gca, 'FontSize',18)
xlabel('\rho'); ylabel('V')

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