设计出完成以下功能的 Matlab 程序：在已知光纤参数和光波长条件下，求解出满 足光纤近似本征方程的 β （相移常数）值，绘制模式电场在光纤横截面上的分布 图。

在MATLAB中设计这样的程序需要使用数值方法和信号处理工具箱，特别是`pdepe`函数用于解偏微分方程。这里是一个基本的步骤和代码示例：

首先，假设我们有一个简化的近似的纤芯模式传播方程，如FDTD模型的简化版本。对于单模光纤，我们可以用下面的数学表达式表示近似本征方程：

\[ \frac{\partial^2 E}{\partial z^2} + k_0^2 n^2(\lambda) \left( \frac{\partial^2 E}{\partial r^2} + \frac{1}{r} \frac{\partial E}{\partial r} - \frac{j}{k_0 r} \beta E \right) = 0 \]

其中 \(E\) 是模式电场，\(z\) 是纤维方向，\(r\) 是径向坐标，\(n(\lambda)\) 是光在光纤中的折射率，\(k_0 = \frac{2\pi}{\lambda}\) 是波数，\(j\) 是虚数单位，\(\beta\) 是我们需要求解的相移常数。

function [beta, Ez] = fiber_solution(n, lambda, L, R)
% n: 光线折射率
% lambda: 光波长
% L: 模式长度
% R: 纤芯半径

% 定义边界条件 (假设边缘全反射)
bc_left = @(r,E) [E(r,0) - E(r,R)];
bc_right = @(r,E) [E(r,0) + E(r,R)];

% 解PDE
[~, Ez, ~] = pdepe('sz', @pdefun, @icfun, [0 R], L, bc_left, bc_right);

% 近似求解β (通常通过数值优化方法)
options = optimoptions('fmincon');
initial_guess = 0; % 初始猜测相移常数
beta = fmincon(@(x) residual(x, n, lambda, Ez), initial_guess, [], [], [], [], options);

% 计算残差函数
residual = @(beta, n, lambda, Ez) norm(diff(Ez, 2, 1));

% 返回β和最后的模式电场分布
end

function dzdx = pdefun(x, t, E, n, lambda, beta)
dzdx = -n^2(lambda) * (diff(E, 2, 1) + diff(E, 1, 2)/x);
dzdx(1,:) = -1i*beta/n(x)*E(1,:);
end

function ic = icfun(x, n, lambda, beta)
% 初始化模式电场，这里可以设置简单的圆柱波或其它形式的初始分布
ic = exp(-100*(x-R).^2); % 示例：高斯分布

运行这个函数后，`beta`将给出满足方程的近似相位常数，而`Ez`则会返回沿光纤长度分布的模式电场。为了可视化横截面的分布，你可以使用`surf`或`pcolor`等函数。