dx=.02; 
xArr=(0:dx:0.5); 
yArr=(0:dx:2); 

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

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

    vMesh(yMesh==0)=(sin(xMesh(yMesh==0)*2*pi/a0));
    vMesh(yMesh==10)=0;

    % iteration step
    for ix=2:length(yArr)-1
        for iy=2:length(xArr)-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')
        colormap hot
        set(gca,'FontSize',18)
        title(['iteration:',num2str(ii)])
        xlabel('x'); ylabel('y'); zlabel('V')
        colorbar
        view(2)
%         caxis([-7 0])
        drawnow 
    end 
end 
%% analytical solution
vArr=vMesh(yMesh==0); 
vAnal=0*vMesh;

for iI=1:20
    sn=sin(xArr*iI*pi/a0); 
    f0=sum(vArr.'.*sn)*dx*2/a0; 
    vAnal=vAnal+f0*sin(xMesh*iI*pi/a0).*exp(-yMesh*iI*pi/a0); 
end 

figure(2)
surf(xMesh,yMesh,(vAnal-vMesh), 'EdgeColor','none')
set(gca,'FontSize',18)
xlabel('x'); ylabel('y'); zlabel('V')
view(2)
colorbar
% caxis([-7 0])

figure(3)
plot(xArr,vArr)
% save('tmp.mat','xMesh','yMesh','vMesh')