w0=180e-6; R=0.01; lmbd=638e-9; k=2*pi/lmbd;
d=.2e-3;
x=-0.0005:0.000001:0.0005;
deltfi=2*pi*d*x/(lmbd*R);
E1=exp(-x.^2/w0.^2);
E2=exp(-(x-d).^2/w0.^2);
figure(gcf);
plot(x,E1,'-r',x,E2,'-g'); pause
P=E1.*E1+E2.*E2+2*cos(deltfi).*E1.*E2;
plot(x,P,'-b');