clear;clc
w0=1;
zp=1;
A0=1;
kw=1;
omg=1;
gamma1=1;
p0=1;
x=(-10:0.01:10)*w0;
z=0:0.01:10;
y=5;%zp*0.5*pi;
[x,z]=meshgrid(x,z);
a=A0*kw/zp./(kw./zp.*cos(z/zp)+2*1i.*sin(z/zp));
b=exp(omg^2*w0*w0.*tan(z/zp)./(tan(z/zp)-1i*kw./zp));
c=cosh(kw./zp*omg*x./(kw./zp.*(cos(z/zp))+2*1i.*sin(z./zp)));
d=cosh(kw./zp*omg*y./(kw./zp.*(cos(z/zp))+2*1i*sin(z./zp)));
e=exp((kw^2*gamma1^2*p0.*tan(z/zp)-2*1i*kw./zp).*(x.^2+y.^2)...
    /w0*w0./(2*1i*kw./zp-4*tan(z/zp)));
phi=a.*b.*c.*d.*e;
mesh(x,z,abs(phi.^2))
view(2)%shading interp