Hi, everyone, why is the laser output power not shown in the figure?

function Pout = dualwavetest4(params,initialConditions, dt, dz, T, L,numSteps)


params = [2825 * 1e-9, 982 * 1e-9, 1620 * 1e-9, 2.136e-21*1e-4, 1.688e-21*1e-4,...
    0.4e-21*1e-4,0.21e-21*1e-4,0.8e-3,6.9e-3,10e-6,120e-6,570e-6,0.37,...
    0.63,0.99,0.00,0.001,0.85,0.006,0.004,0.14,0.34,0.012,0.015,0.18,...
    0.44,1.0,2,2,0.1,0.16,2.8918e-10,1.28,0.362,0.542,1.55,1.12e-21*1e-4,...
    2.31e-21*1e-4,3.02e-21*1e-4,3.2e+026,0.03,0.01625,0.0325,0.7511,0.0024,...
    0.052,0,0.04,0,0.04,0.994,0.04,16,3e8,6.626e-34,10,5];

%params(1) = lambda_s = 2825 * 1e-9; 
%params(2) = lambda_p1 = 982 * 1e-9; 
%params(3) = lambda_p2 = 1620 * 1e-9; 
%params(4) = sigma_ap1 = 2.136e-21*1e-4; 
%params(5) = sigma_ep1 = 1.688e-21*1e-4; 
%params(6) = sigma_ap2 = 0.4e-21*1e-4;
%params(7) = sigma_ep2 = 0.21e-21*1e-4; 
%params(8) = tau1 = 0.8e-3; 
%params(9) = tau2 = 6.9e-3;
%params(10)= tau3 = 10e-6;
%params(11)= tau4 = 120e-6;
%params(12)= tau5 = 570e-6;
%params(13)= beta21 = 0.37;
%params(14)= beta20 = 0.63;
%params(15)= beta32 = 0.99;
%params(16)= beta31 = 0.00;
%params(17)= beta30 = 0.001;
%params(18)= beta43 = 0.85;
%params(19)= beta42 = 0.006;
%params(20)= beta41 = 0.004;
%params(21)= beta40 = 0.14;
%params(22)= beta54 = 0.34;
%params(23)= beta53 = 0.012;
%params(24)= beta52 = 0.015;
%params(25)= beta51 = 0.18;
%params(26)= beta50 = 0.44;
%params(27)= beta10 = 1.0;
%params(28)= g1 = 2;
%params(29)= g2 = 2;
%params(30)= b1 = 0.1;
%params(31)= b2 = 0.16;
%params(32)= Aeff = 2.8918e-10;
%params(33)= W11 = 1.28;
%params(34)= W22 = 0.362;
%params(35)= W50 = 0.542;
%params(36)= W42 = 1.55;
%params(37)= sigma_25 = 1.12e-21*1e-4; 
%params(38)= sigma_52 = 2.31e-21*1e-4;

%params(39)= sigma_es = 3.02e-21*1e-4;% 
%     
%     
%params(40)= NEr = 3.2e+026;
%params(41)= alpha_p1 = 0.03;
%params(42)= alpha_p2 = 0.01625;
%params(43)= alpha_s = 0.0325;
%params(44)= gamma_s = 0.7511;
%params(45)= gamma_p1 = 0.0024; % 
%params(46)= gamma_p2 = 0.052; % 
%params(47)= Rp1 = 0;
%params(48)= Rp2 = 0.04;
%params(49)= Rp3 = 0;
%params(50)= Rp4 = 0.04;
%params(51)= Rs1 = 0.994;
%params(52)= Rs2 = 0.04;
%params(53)= L = 16;

%params(54)= c = 3e8;
%params(55)= h = 6.626e-34;

%Ppl = [Pinput_p1; Pinput_p2];  
%params(56)= Pinput_p1 = 10;  
%params(57)= Pinput_p2 = 5;  
%params(58)= Ppr = 0;   


N = zeros(6, 1); % 
y = [params(56);0;params(57);0;0;0]; % 
Z = linspace(0,params(53),100); % 
t = linspace(0, 90, 100); % 
numSteps = length(t);
dt = 0.01;
dz = 0.01;


y_history = zeros(length(y), numSteps); 
N_history = zeros(6, numSteps); 

for i = 1:numSteps
    
    for j = 1:length(Z)-1%   
        
        %             y = apply_boundary_conditions(y, params, Z(j));
        
        
        %   fprintf('Boundary conditions at Z=%.2f applied. Before: Pump powers: %.2e, %.2e; Laser power: %.2e\n', ...
        %        Z,y(1), y(3), y(5));
        
        
        [N, y] = rk4_step(params, N, y, Z(j), dt, dz);
    end
    
    
    y_history(:, i) = y;
    N_history(:, i) = N;
    
    
    
    
    %         y = apply_boundary_conditions(y, params, 0); 
    %         y = apply_boundary_conditions(y, params, params(53));  
    
    %     fprintf('After: Pump powers: %.2e, %.2e; Laser power: %.2e\n', ...
    %        y_history(1,:), y_history(3,:), y_history(5,:));
    
end


figure;
subplot(2,1,1);
plot(Z, y_history(1,:), 'b.-', 'DisplayName', 'Pp1+');
hold on;
plot(Z, y_history(2, :), 'g*-', 'DisplayName', 'Pp1-');
plot(Z, y_history(3, :), 'p-', 'DisplayName', 'Pp2+');
plot(Z, y_history(4, :), 'y>-', 'DisplayName', 'Pp2-');
plot(Z, y_history(5, :), 'r', 'DisplayName', 'Ps+');
plot(Z, y_history(6, :), 'k--', 'DisplayName', 'Ps-');
grid on;
title('Pump and laser powers');
legend('Location', 'best');
xlabel('Position z (m)');
ylabel('Power (W)');


subplot(2,1,2);
for k = 1:6
    plot(Z, N_history(k, :), 'DisplayName', strcat('energy level ', num2str(k)));
    hold on;
end
grid on;
title('Relative population density');
legend('Location', 'best');
xlabel('Position z (m)');
ylabel('N/N');






end

function y = apply_boundary_conditions(y, params, Z)


%     if Z == 0 || Z == params(53) 
%         % Z = 0 
%         y(1) = params(47) * y(2) + params(56); % Pp1
%         y(3) = params(49) * y(4) + params(57); % Pp2
%         
%         %  Z = L
%         y(2) = params(48) * y(1); 
%         y(4) = params(50) * y(3); 
%         
%         % Boundary conditions
%         y(5) = params(51) * y(6); % Emission conditions of laser power
%         y(6) = params(52) * y(5); 
%     end


end

function [N, y] = rk4_step(params, N, y, Z, dt, dz)
params = [2825 * 1e-9, 982 * 1e-9, 1620 * 1e-9, 2.136e-21*1e-4, 1.688e-21*1e-4,...
    0.4e-21*1e-4,0.21e-21*1e-4,0.8e-3,6.9e-3,10e-6,120e-6,570e-6,0.37,...
    0.63,0.99,0.00,0.001,0.85,0.006,0.004,0.14,0.34,0.012,0.015,0.18,...
    0.44,1.0,2,2,0.1,0.16,2.8918e-10,1.28,0.362,0.542,1.55,1.12e-21*1e-4,...
    2.31e-21*1e-4,3.02e-21*1e-4,3.2e+026,0.03,0.01625,0.0325,0.7511,0.0024,...
    0.052,0,0.04,0,0.04,0.994,0.04,16,3e8,6.626e-34,10,5];

k1_N = dt * rate_equations(N, y, params);
k1_y = dz * transmission_equations(N, y, params, Z);

k2_N = dt * rate_equations(N + 0.5 * k1_N, y + 0.5 * k1_y, params);
k2_y = dz * transmission_equations(N + 0.5 * k1_N, y + 0.5 * k1_y, params, Z);

k3_N = dt * rate_equations(N + 0.5 * k2_N, y + 0.5 * k2_y, params);
k3_y = dz * transmission_equations(N + 0.5 * k2_N, y + 0.5 * k2_y, params, Z);

k4_N = dt * rate_equations(N + k3_N, y + k3_y, params);
k4_y = dz * transmission_equations(N + k3_N, y + k3_y, params, Z);

N = N + (1/6) * (k1_N + 2*k2_N + 2*k3_N + k4_N);
y = y + (1/6) * (k1_y + 2*k2_y + 2*k3_y + k4_y);


end


function dNdt = rate_equations(N, y, params)
params = [2825 * 1e-9, 982 * 1e-9, 1620 * 1e-9, 2.136e-21*1e-4, 1.688e-21*1e-4,...
    0.4e-21*1e-4,0.21e-21*1e-4,0.8e-3,6.9e-3,10e-6,120e-6,570e-6,0.37,...
    0.63,0.99,0.00,0.001,0.85,0.006,0.004,0.14,0.34,0.012,0.015,0.18,...
    0.44,1.0,2,2,0.1,0.16,2.8918e-10,1.28,0.362,0.542,1.55,1.12e-21*1e-4,...
    2.31e-21*1e-4,3.02e-21*1e-4,3.2e+026,0.03,0.01625,0.0325,0.7511,0.0024,...
    0.052,0,0.04,0,0.04,0.994,0.04,16,3e8,6.626e-34,10,5];


%   RGSA = params(2)*params(45)*(params(4)*N(0)-params(5)*N(2))*(y(1)+y(2))/(params(55)*params(54)*params(32));
%   RESA1 = params(2)*params(45)*(params(37)*N(2)-params(38)*N(5))*(y(1)+y(2))/(params(55)*params(54)*params(32));
%   RESA2 = params(3)*params(46)*(params(6)*N(1)-params(7)*N(3))*(y(3)+y(4))/(params(55)*params(54)*params(32));
%   RSE = params(1)*params(44)*params(39)*(b2*N(2)-g2/g1*b1*N(1))*(y(5)+y(6))/(params(55)*params(54)*params(32));


dNdt = zeros(6,1);
dNdt(6) = -N(6)/params(6)+params(6)*N(3)*N(3)-params(35)*N(6)*N(1)+params(36)*N(5)*N(3)+params(2)*params(45)*(params(37)*N(3)-params(38)*N(6))*(y(1)+y(2))/(params(55)*params(54)*params(32));
dNdt(5) = -N(5)/params(11)-params(36)*N(5)*N(3)+params(22)*N(6)/params(12);
dNdt(4) = -N(4)/params(10)+params(33)*N(2)*N(2)+params(35)*N(6)*N(1)+params(18)/params(11)*N(5)+params(23)/params(12)*N(6)+params(3)*params(46)*(params(6)*N(2)-params(7)*N(4))*(y(3)+y(4))/(params(55)*params(54)*params(32));
dNdt(3) = -N(3)/params(9)-2*params(34)*N(3)*N(3)-params(36)*N(5)*N(3)+params(15)/params(10)*N(4)+params(19)/params(11)*N(5)+params(24)/params(12)*N(6)+params(2)*params(45)*(params(4)*N(1)-params(5)*N(3))*(y(1)+y(2))/(params(55)*params(54)*params(32))-(params(2)*params(45)*(params(37)*N(3)-params(38)*N(6))*(y(1)+y(2))/(params(55)*params(54)*params(32)))-(params(1)*params(44)*params(39)*(params(31)*N(3)-params(29)/params(28)*params(30)*N(2))*(y(5)+y(6))/(params(55)*params(54)*params(32)));
dNdt(2) = -N(2)/params(8)-2*params(33)*N(2)*N(2)+params(35)*N(6)*N(1)+params(36)*N(5)*N(3)+params(13)/params(9)*N(3)+params(16)/params(10)*N(4)+params(20)/params(11)*N(5)+params(25)/params(12)*N(6)+params(1)*params(44)*params(39)*(params(31)*N(3)-params(29)/params(28)*params(30)*N(2))*(y(5)+y(6))/(params(55)*params(54)*params(32))-(params(3)*params(46)*(params(6)*N(2)-params(7)*N(4))*(y(3)+y(4))/(params(55)*params(54)*params(32)));
dNdt(1) = params(33)*N(2)*N(2)+params(34)*N(3)*N(3)-params(35)*N(6)*N(1)+params(27)/params(8)*N(2)+params(14)/params(9)*N(4)+params(17)/params(10)*N(3)+params(21)/params(11)*N(5)+params(26)/params(12)*N(6)-(params(2)*params(45)*(params(4)*N(1)-params(5)*N(3))*(y(1)+y(2))/(params(55)*params(54)*params(32)));
%N(6)+N(1)+N(2)+N(3)+N(4)+N(5)= params(40);




end


function dydz = transmission_equations(N, y, params, Z)
params = [2825 * 1e-9, 982 * 1e-9, 1620 * 1e-9, 2.136e-21*1e-4, 1.688e-21*1e-4,...
    0.4e-21*1e-4,0.21e-21*1e-4,0.8e-3,6.9e-3,10e-6,120e-6,570e-6,0.37,...
    0.63,0.99,0.00,0.001,0.85,0.006,0.004,0.14,0.34,0.012,0.015,0.18,...
    0.44,1.0,2,2,0.1,0.16,2.8918e-10,1.28,0.362,0.542,1.55,1.12e-21*1e-4,...
    2.31e-21*1e-4,3.02e-21*1e-4,3.2e+026,0.03,0.01625,0.0325,0.7511,0.0024,...
    0.052,0,0.04,0,0.04,0.994,0.04,16,3e8,6.626e-34,10,5];


dydz = zeros(6,1);
dydz(1) = (-params(45)*(params(4)*N(1)-params(5)*N(3)+params(37)*N(3)-params(38)*N(6))-params(41))*y(1);
dydz(2) = (params(45)*(params(4)*N(1)-params(5)*N(3)+params(37)*N(3)-params(38)*N(6))+params(41))*y(2);
dydz(3) = (-params(46)*(params(6)*N(2)-params(7)*N(4))-params(42))*y(3);
dydz(4) = (params(46)*(params(6)*N(2)-params(7)*N(4))+params(42))*y(4);
dydz(5) = (params(44)*params(39)*(params(31)*N(3)-params(29)/params(28)*params(30)*N(2))-params(43))*y(5);
dydz(6) = (-params(44)*params(39)*(params(31)*N(3)-params(29)/params(28)*params(30)*N(2))+params(43))*y(6);




fprintf('At Z=%.2f, Rate of change for pump powers: %.2e, %.2e; Laser power: %.2e\n', ...
    Z, dydz(1), dydz(3), dydz(5));


if Z == 0
    %          Z = 0 
    y(1) = params(47)*y(2)+params(56);%y(1) = Rp1*y(2)+P_input_p1;  
    y(3) = params(49)*y(4)+params(57);%y(3) = Rp3*y(4)+P_input_p2;
    y(5) = params(51)*y(6);%y(5) = Rs1*y(6);
    
    
elseif Z == params(53)
    y(2) = params(48)*y(1);%y(2) = Rp2*y(1);         
    y(4) = params(50)*y(3);%y(4) = Rp4*y(3);
    y(6) = params(52)*y(5);%y(6) = Rs2*y(5);
    
    
end
fprintf('At Z=%.2f, Rate of change for pump powers: %.2e, %.2e; Laser power: %.2e\n', ...
    Z, dydz(1), dydz(3), dydz(5));

end

[dpb]

At least format the code first...

Answer 1

[Walter Roberson]

plot(Z, y_history(4, :), 'y>-', 'DisplayName', 'Pp2-');
plot(Z, y_history(5, :), 'r', 'DisplayName', 'Ps+');
plot(Z, y_history(6, :), 'k--', 'DisplayName', 'Ps-');

All of those values are identical to zero.

dNdt(6) = -N(6)/params(6)+params(6)*N(3)*N(3)-params(35)*N(6)*N(1)+params(36)*N(5)*N(3)+params(2)*params(45)*(params(37)*N(3)-params(38)*N(6))*(y(1)+y(2))/(params(55)*params(54)*params(32));

Every term there depends upon a value from N. When N is initialized to all zero, the result must be identical to zero.

N = N + (1/6) * (k1_N + 2*k2_N + 2*k3_N + k4_N);

all-zero plus terms that are all-zero, gives all zero. So N stays all zero.

dydz(4) = (params(46)*(params(6)*N(2)-params(7)*N(4))+params(42))*y(4);
dydz(5) = (params(44)*params(39)*(params(31)*N(3)-params(29)/params(28)*params(30)*N(2))-params(43))*y(5);
dydz(6) = (-params(44)*params(39)*(params(31)*N(3)-params(29)/params(28)*params(30)*N(2))+params(43))*y(6);

y(4), y(5), y(6) happen to all be 0, and due to the multiplication by y(4), y(5), y(6), those derivatives all come out as 0. So the last 3 y values all come out as 0.

[LinHai]

Thanks for your answer. Here, I set the initial value of the laser power to 2 W, logically speaking, the laser power will increase and stabilize to a fixed value, why the laser power keeps decreasing? And how should I set the initial value of N?