x=0.00001:.01:.5;  
y=0.00001:.01:.5;  
[xx,yy]=meshgrid(x,y);  

a0 = 1;  
q = 1;  
a = 2/a0;  

coord = (xx-0).^2+(yy-0).^2; 
r = sqrt(coord);  

phi = q * exp((-a).*r) ./ r .*(1.+a.*r./2);  
phi1 = q * exp((-a).*r) ./ r;  
phi2 = q * exp((-a).*r) ./ r .* a .* r./ 2;  

lp = del2(phi); 
mesh(lp);  
title('Graph of the Laplacian of Phi') 
xlabel('xdistance') 
ylabel('y distance') 
zlabel('laplacian of Phi') 

rho = lp ./ ((-4) .* pi)

mesh(rho) 
title('The charge distribution') 
xlabel('x distance') 
ylabel('y distance')
zlabel('Rho')