function Pout = fiberlaser_twoend
global R1 R2 Ppl Ppr sigma_ap sigma_ep sigma_as sigma_es gamma_s ...
gamma_p N alpha_p alpha_s Pssat Ppsat Ppsp Ppsv mu kappa elta
%参数设置
lambda_s = 1530 * 1e-9;%发射波长m
lambda_p = 978 * 1e-9;%泵浦波长m
tau = 10e-3;%上能及寿命Er离子亚能级10ms
sigma_ap = 7.9e-24;%泵浦光的吸收截面26e-21*1e-4
sigma_ep = 6.5e-24;%泵浦光的发射截面26e-21*1e-4
sigma_as = 3e-27;%光纤激光的吸收截面
sigma_es = 2.5e-25;%光纤激光的发射截面
A_c = 28.e-12;%纤芯的截面积对应6um
N = 2.5e25;%掺杂粒子数密度700ppm
alpha_p = 2e-3;%光纤对泵浦光的损耗
alpha_s = 4e-4;%光纤对激光的损耗
gamma_s = 0.82;%激光功率填充因子
gamma_p = 0.0024;%泵浦光功率填充因子
R1 =0;%前反射镜反射率
R2 =1;%后反射镜反射率(出光方向)
L = 10;%光纤长度
%物理常数及中间过程参数计算
c = 3e8;%光速
h = 6.626e-34;%普朗克常数
nu_s = c/lambda_s;%发射光频率
nu_p = c/lambda_p;%泵浦光频率
Pssat = h * nu_s * A_c/( gamma_s * (sigma_es+sigma_as) * tau);
Ppsat = h * nu_p * A_c/( gamma_p * (sigma_ep+sigma_ap) * tau);
%抽运光功率设置
Ppl = 50e-3;
Ppr = 0.0;
%`端面抽运的光纤激光器边值问题数值求解
OPTION = bvpset('Stats','ON');
solinit = bvpinit(linspace(0,L,10),[Ppl Ppr 30 Ppr]);
sol = bvp4c(@f,@fsbc,solinit);
%数值计算结果分析和显示
x = [sol.x];
y = [sol.y];
nz = [(sigma_ap/(sigma_ap+sigma_ep)*(y(1,:)+y(2,:))/Ppsat+...
sigma_as/(sigma_as+sigma_es)*(y(3,:)+y(4,:))/Pssat)./...
((y(1,:)+y(2,:))/Ppsat+1+(y(4,:)+y(3,:))/Pssat)];
Pout = y(3,end)*(1-R2);
figure;
subplot(2,1,1,'LineWidth', 1.5)
plot(x,y(1,:)*1000,'m.-',x,y(2,:)*1000,'b.-',x,y(3,:)*1000,'r',x,y(4,:)*1000,'k--','LineWidth', 1.5);
grid on;
title('Pump and laser powers', 'FontName', 'Times New Roman', 'FontSize', 15);
legend('Pp+(z)','Pp-(z)','Ps+(z)','Ps-(z)', 'FontName', 'Times New Roman', 'FontSize', 10);
xlabel('Position z (m)', 'FontName', 'Times New Roman', 'FontSize', 15);
ylabel('Power (mW)', 'FontName', 'Times New Roman', 'FontSize', 15);
subplot(2,1,2,'LineWidth', 1.5)
plot(x,nz,'LineWidth', 1.5)
grid on;
title('Relative population density', 'FontName', 'Times New Roman', 'FontSize', 15)
xlabel('Position z (m)', 'FontName', 'Times New Roman', 'FontSize', 15);
ylabel('N_2/N', 'FontName', 'Times New Roman', 'FontSize', 15);
%端面抽运的光纤激光器速率方程组
function dy = f(x,y)
global sigma_ap sigma_ep sigma_as sigma_es gamma_s gamma_p N ...
alpha_p alpha_s Pssat Ppsat Ppsp Ppsv
dy = zeros(4,1);
ep=1e-5;%瑞利散射分布反馈
es=4e-5;
N21=N*(sigma_ap/(sigma_ap+sigma_ep)*(y(1)+y(2))/Ppsat+sigma_as/...
(sigma_as+sigma_es)*(y(3)+y(4))/Pssat)/((y(1)+y(2))/Ppsat+1+(y(4)+y(3))/Pssat);
dy(1)=(-gamma_p*(sigma_ap*N-(sigma_ap+sigma_ep)*N21)-alpha_p)*y(1)+ep*y(2);
dy(2)=-(-gamma_p*(sigma_ap*N-(sigma_ap+sigma_ep)*N21)-alpha_p)*y(2)-ep*y(1);
dy(3)=(gamma_s*((sigma_as+sigma_es)*N21-sigma_as*N)-alpha_s)*y(3)+es*y(4);
dy(4)=-(gamma_s*((sigma_as+sigma_es)*N21-sigma_as*N)-alpha_s)*y(4)-es*y(3);
%端面抽运的光纤激光器边界条件
function res = fsbc(y0,yL)
global R1 R2 Ppl Ppr
res = [y0(1)-Ppl
yL(2)-Ppr
y0(3)-R1*y0(4)
yL(4)-R2*yL(3)];
怎么把上面MATLAB代码画出的曲线导出在Excel里面