拉曼放大器 matlab,求助各位老师，ode45计算拉曼放大器结果不收敛

function Pout = fiberlaser_raman

global lambda_0 lambda_1 lambda_2 lambda_3 D_core D_clad1 D_clad2 D_clad3 D_clad4 S_core S_clad4 S_clad3 S_clad2 S_clad1 g_Rcore g_R1 g_R2 g_R3 g_R4 ...

g_1 g_2 g_3 alpha_0 alpha_1 alpha_2 alpha_3

lambda_0 = 976 * 1e-9;

lambda_1 = 1070 * 1e-9;

lambda_2 = 1181 * 1e-9;

lambda_3 = 1318 * 1e-9;

D_core = 5 * 1e-6;% 1318

D_clad4 = 5.5 * 1e-6;% 1181

D_clad3 = 14 * 1e-6;% 1070

D_clad2 = 27 * 1e-6;%此层为泵光

D_clad1 = 52.5 * 1e-6;%外束缚用的包层

S_core = (D_core/2)^2*pi;

S_clad4 = (D_clad4/2)^2*pi;

S_clad3 = (D_clad3/2)^2*pi;

S_clad2 = (D_clad2/2)^2*pi;

S_clad1 = (D_clad1/2)^2*pi;

g_Rcore =  5.3 * 1e-14;

g_R1 = 5.3 * 1e-14;

g_R2 = 5.3 * 1e-14;

g_R3 = 5.3 * 1e-14; %后期再改系数，与Ge掺杂浓度相关

g_R4 = 5.3 * 1e-14;

g_1 = g_Rcore*S_core/S_clad3+g_R4*(S_clad4-S_core)/S_clad3+g_R3*(S_clad3-S_clad4)/S_clad3;

g_2 = g_Rcore*S_core/S_clad4+g_R4*(S_clad4-S_core)/S_clad4;

g_3 = g_R3;

% g_1 = g_R1;

% g_2 = g_R2;

% g_3 = g_R3;

alpha_0 = 0.0005;

alpha_1 = 0.0003;

alpha_2 = 0.0003;

alpha_3 = 0.0003;

L = 100;

Ppl = 2000;

Ppr1 = 50;

Ppr2 = 0.1;

Ppr3 = 0.01;

options = odeset('RelTol',1e-5,'AbsTol',[1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5]);

[x,y] = ode45(@fiberlaser,[0 L],[Ppl 0 Ppr1 0 Ppr2 0 Ppr3],options);

%%计算结果显示与分析

% x = [sol.x];

% y = [sol.y];

Pout = y(7,end);

figure;

plot(x,y(:,1),'b.-',x,y(:,3),'g*-',x,y(:,5),'r',x,y(:,7),'k--');

grid on;

title('Pump and raman output');

legend('Pp','P1st','P2nd','P3rd')

xlabel('Position z (m)');

ylabel('Power (w)');

% 拉曼转换方程组

function dy = fiberlaser(x,y)% x 为光纤长度，y 为 不同包层光功率

global lambda_0 lambda_1 lambda_2 lambda_3 S_core S_clad4 S_clad3 S_clad2 S_clad1 g_Rcore g_R1 g_R2 g_R3 g_1 g_2 g_3 ...

alpha_0 alpha_1 alpha_2 alpha_3

dy = zeros(7,1);

dy(1) = -lambda_1*g_1/(lambda_0*S_clad1)*[y(3)*y(1)+y(2)*y(1)]-alpha_0*y(1);

dy(2) = g_1*y(1)*y(2)/S_clad1-lambda_2*y(2)*g_2/(lambda_1*S_clad1)*[y(4)+y(5)]-alpha_1*y(2);

dy(3) = g_1*y(1)*y(3)/S_clad1-lambda_2*y(3)*g_2/(lambda_1*S_clad1)*[y(4)+y(5)]-alpha_1*y(3);

dy(4) = g_2*y(2)*y(4)/S_clad1+g_2*y(3)*y(4)/S_clad2-lambda_3*g_3/(lambda_2*S_clad2)*y(4)*(y(6)+y(7))-alpha_2*y(4);

dy(5) = g_2*y(2)*y(5)/S_clad1+g_2*y(3)*y(5)/S_clad2-lambda_3*g_3/(lambda_2*S_clad3)*y(5)*(y(6)+y(7))-alpha_2*y(5);

dy(6) = g_3*y(4)*y(6)/S_clad2-g_3*y(5)*y(6)/S_clad3-alpha_3*y(6);

dy(7) = g_3*y(4)*y(7)/S_clad2-g_3*y(5)*y(7)/S_clad3-alpha_3*y(7);