Mathematica 中将 高 斯 波 ：u[x, 0] == Exp[-x^2/(b*w^2)]作 初 值 时 ， 求 解 衍 射 傍 轴 方 程 ：I*D[u[x, z], z] + (1/(2*2))*D[u[x, z], {x, 2}] == 0（k=2*Pi/\[Lambda]）。比较 b=1,b=2 和 b=4 时强度随着传输距离 z 增加 演化图，距离要求 z>k*b*w^2

以下是 Mathematica 的代码实现：

λ = 1;  (*波长*)
k = 2*π/λ;  (*波数*)
p = 10;  (*计算区间*)
steps = 100;  (*步数*)
dz = p/steps;  (*步长*)

bValues = {1, 2, 4};  (*高斯波参数*)
w = 1;  (*高斯波参数*)

sols = Table[
   NDSolve[{I*D[u[x, z], z] + (1/(2*b^2))*D[u[x, z], {x, 2}] == 0, u[x, 0] == Exp[-x^2/(b*w^2)]}, u, {x, -p, p}, {z, 0, p}], 
   {b, bValues}
];

frames = Table[
   Plot[Abs[u[x, z] /. sols[[i]]]^2, {x, -p, p}, PlotRange -> {0, 1}, PlotLabel -> StringJoin["b=", ToString[bValues[[i]]]]],
   {i, Length[bValues]},
   {z, 0, p, dz}
];

ListAnimate[frames, AnimationRate -> 5]

其中，`Table` 函数用于循环求解不同 b 值下的衍射傍轴方程，`Plot` 函数用于绘制结果，`ListAnimate` 函数用于生成动画。最终，我们得到了不同 b 值下衍射傍轴方程的演化图。需要注意的是，由于数值计算的限制，演化图只能在距离 z > k*b*w^2 时才有意义。