以矩形腔为例，用matlab说明 Fox-Li 数值迭代法计算时 的具体形式。

Fox-Li 数值迭代法是用于计算光学腔模式的一种常用方法。下面以矩形腔为例，介绍 Fox-Li 数值迭代法的具体形式。

假设矩形腔的两个短边长度为 $a$，两个长边长度为 $b$。设矩形腔内的电磁场分别为 $E_x$、$E_y$、$H_x$ 和 $H_y$，则它们满足以下波动方程：

$$
\nabla^2 E_x+k_0^2n^2E_x=0 \tag{1}
$$

$$
\nabla^2 E_y+k_0^2n^2E_y=0 \tag{2}
$$

$$
\nabla^2 H_x+k_0^2n^2H_x=0 \tag{3}
$$

$$
\nabla^2 H_y+k_0^2n^2H_y=0 \tag{4}
$$

其中 $k_0=\frac{2\pi}{\lambda}$，$n$ 为介质折射率，$\lambda$ 为光波长。

将电场和磁场表示为横向和纵向的模式，即 $E_x=E_{x0}\cos(k_xx)\cos(k_yy)$，$E_y=E_{y0}\sin(k_xx)\sin(k_yy)$，$H_x=H_{x0}\sin(k_xx)\sin(k_yy)$，$H_y=H_{y0}\cos(k_xx)\cos(k_yy)$，代入波动方程中得到：

$$
\frac{\partial^2E_{x0}}{\partial x^2}+\frac{\partial^2E_{x0}}{\partial y^2}+k_0^2n^2E_{x0}=0 \tag{5}
$$

$$
\frac{\partial^2E_{y0}}{\partial x^2}+\frac{\partial^2E_{y0}}{\partial y^2}+k_0^2n^2E_{y0}=0 \tag{6}
$$

$$
\frac{\partial^2H_{x0}}{\partial x^2}+\frac{\partial^2H_{x0}}{\partial y^2}+k_0^2n^2H_{x0}=0 \tag{7}
$$

$$
\frac{\partial^2H_{y0}}{\partial x^2}+\frac{\partial^2H_{y0}}{\partial y^2}+k_0^2n^2H_{y0}=0 \tag{8}
$$

利用 Fox-Li 方法，我们可以将矩形腔内的电磁场表示为以下形式：

$$
E_{x0}(x,y)=\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}a_{mn}\sin(\frac{m\pi x}{a})\cos(\frac{n\pi y}{b}) \tag{9}
$$

$$
E_{y0}(x,y)=\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}b_{mn}\cos(\frac{m\pi x}{a})\sin(\frac{n\pi y}{b}) \tag{10}
$$

$$
H_{x0}(x,y)=\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}c_{mn}\cos(\frac{m\pi x}{a})\sin(\frac{n\pi y}{b}) \tag{11}
$$

$$
H_{y0}(x,y)=\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}d_{mn}\sin(\frac{m\pi x}{a})\cos(\frac{n\pi y}{b}) \tag{12}
$$

将式（9）~（12）代入式（5）~（8）中，利用正交性和归一化条件，可以得到以下的本征值方程组：

$$
\frac{\omega_{mn}^2}{c^2}a_{mn}=-\frac{\pi^2}{ab}(m^2+n^2)a_{mn}+\frac{\pi^2}{ab}\sum_{p=1}^{\infty}\sum_{q=1}^{\infty}\epsilon_{mp}\epsilon_{nq}(\frac{m}{a}\frac{p}{a}+\frac{n}{b}\frac{q}{b})b_{pq} \tag{13}
$$

$$
\frac{\omega_{mn}^2}{c^2}b_{mn}=-\frac{\pi^2}{ab}(m^2+n^2)b_{mn}+\frac{\pi^2}{ab}\sum_{p=1}^{\infty}\sum_{q=1}^{\infty}\epsilon_{mp}\epsilon_{nq}(\frac{m}{a}\frac{p}{a}+\frac{n}{b}\frac{q}{b})a_{pq} \tag{14}
$$

$$
\frac{\omega_{mn}^2}{c^2}c_{mn}=-\frac{\pi^2}{ab}(m^2+n^2)c_{mn}+\frac{\pi^2}{ab}\sum_{p=1}^{\infty}\sum_{q=1}^{\infty}\epsilon_{mp}\epsilon_{nq}(\frac{m}{a}\frac{p}{a}+\frac{n}{b}\frac{q}{b})d_{pq} \tag{15}
$$

$$
\frac{\omega_{mn}^2}{c^2}d_{mn}=-\frac{\pi^2}{ab}(m^2+n^2)d_{mn}+\frac{\pi^2}{ab}\sum_{p=1}^{\infty}\sum_{q=1}^{\infty}\epsilon_{mp}\epsilon_{nq}(\frac{m}{a}\frac{p}{a}+\frac{n}{b}\frac{q}{b})c_{pq} \tag{16}
$$

其中 $\omega_{mn}$ 为模式频率，$c$ 为光速，$\epsilon_{mn}$ 为 Kronecker delta 函数，当 $m=n$ 时 $\epsilon_{mn}=1$，否则 $\epsilon_{mn}=0$。

本征值方程组可以用 Fox-Li 数值迭代法进行求解。具体来说，我们可以采用以下步骤：

1. 给定初始的 $a_{mn}$、$b_{mn}$、$c_{mn}$ 和 $d_{mn}$ 的估计值。

2. 根据式（13）~（16）计算下一次迭代的 $a_{mn}$、$b_{mn}$、$c_{mn}$ 和 $d_{mn}$ 的值。

3. 判断迭代是否收敛。如果没有收敛，则回到步骤2；如果收敛，则得到模式频率 $\omega_{mn}$ 和对应的模式场分布。

4. 重复步骤1~3，直到求得所有的模式频率和模式场分布。

以上就是 Fox-Li 数值迭代法计算矩形腔模式时的具体形式。在实际计算中，我们通常采用 MATLAB 等数值计算软件进行计算。