[lv19860312]

程序如下：
syms x y z w0 lad p l
x=-0.005:0.0001:0.005;
y=-0.005:0.0001:0.005;
z=-100:2:100;
w0=2*10^-3;
lad=632*10^-9;
Z0=pi*w0^2/lad;
k=2*pi/lad;
w=w0*sqrt(1+(z./Z0));
R=Z0*(z./Z0+Z0./z);
t=2*(x.^2+y.^2)./w.^2;
p=1;l=1;
L=LaguerreL(p,l,t);
u=sqrt(1+(z/Z0))*.L.*(t.^l/2).*exp(-0.5*t).*exp(-i*k.*(z+(x.^2+y.^2)/2*R)...
-i*(2*p+l+1)*(tan(z./Z0)).^-1).*exp(-i*l.*atan(x./y));
p=patch(isosurface(x,y,z,u,0));
daspect([1 1 1])
view(3)
camlight; lighting phong
其中还调用了下面的函数
function L=LaguerreL(p,l,t)
%拉盖尔多项式
maple('LaguerreL',p,l,t);
L=maple('simplify',ans,'LaguerreL');