clear all
Q1=5e-12; 
e0=8.854e-12;
k=1/(4*pi*e0);  
r1=[1, 2]*1e-2;

dx=0.5e-3;dy=0.5e-3; maxx=4e-2;maxy=4e-2;
x=0:dx:maxx; y=0:dy:maxy;
[X,Y]=meshgrid(0:dx:maxx,0:dy:maxy);


for i=1:1:length(x)
   for j=1:1:length(y)
      p1x=r1(1,1)-x(i);
      p1y=r1(1,2)-y(j);
      psqr=sqrt(p1x^2+p1y^2);
      if psqr==0
         psqr=1e-10;
      end
      E1x(i,j)=k*Q1*p1x/psqr^3;
      E1y(i,j)=k*Q1*p1y/psqr^3;
   end
end  
Eabs=(sqrt(E1x.^2+E1y.^2));
Emax=max(max(Eabs));
[i,j]=find(Eabs==0);
Eabs(i,j)=Emax;

figure,mesh(X,Y,Eabs)