Why the plot doesn't show the result? I am trying to Solve Q-Switched Fiber Laser Rate Equation

Hi everyone

I have no problem with error line codes. However my problem is when I try to plot the result, it doesn't show anything, just straight line at y = 0.

This is my code, thanks for your help.

clear
clc
ti = 0; 
tf = 3E-8; 
tspan=[ti tf]; 
y0=[0; 0; 0]; 
%y(1) = Photon Density
%y(2) = Inverted Population Density
%y(3) = Photocarrier Density
[t,y] = ode45(@rate_eq,tspan,y0);
figure(1); 
hold on 
plot(t,y(:,1)); 
title('Photon Density'); 
xlabel('Time'); 
ylabel('Photon Density'); 
hold off 
figure(2); 
hold on 
plot(t,y(:,2)); 
title('Inverted Population Density '); 
xlabel('Time'); 
ylabel('Inverted Population Density'); 
hold off 
figure(3); 
hold on 
plot(t,y(:,3)); 
title('Photocarrier Density'); 
xlabel('Time'); 
ylabel('Photocarrier Density'); 
hold off 
 
function dy = rate_eq( t,y )
%Rate equation for Q-switched fiber laser.
dy = zeros(3,1);
% Planck constant (m^2kg/s)
h = 6.62606957E-34;
% Speed of light
c = 299792458;
% Dissipative optical loss
delta = 0.4;
%Inversion reduction factor
gamma = 1.8;
% Length of total fiber (m)
lr = 15;
% Length of active fiber (m)
L = 3;
% Thickness of the SA (m)
Lsa = 1E-5;
% Output coupling ratio
R = 0.95;
% Stimulated emission area (m2)
sigmaes = 2E-25;
% Spontataneous decay time(t)
tg = 1E-2;
% Saturation photocarrier density 
Nsa = 2.4E27;
% SA carrier recombination time (s) 
tsa = 0.1E-9;
% nonsaturable loss (s) 
lambdans = 0.4;
% modulation depth 
lambdas = 0.1;
% Effective doping area of active fiber (m2) 
A = 1.26E-11;
% Pumping power (W)
Pp = 2E-2;
% Wavelength of signal light (m)
lambdaP = 1530E-9;
% Frequency of signal (Hz)
v = c/lambdaP;
% Round trip transit time (Hz)
tr = lr/c;
% Pumprate of active medium
Wp = Pp/(h*v*A*L);
% Saturable absorption
lambdana = (lambdas/(1+(y(3)/Nsa)))+lambdans;
%lambdana = 0.5;
% Change of Photon density
dy(1) = (y(1)/tr)*(2*sigmaes*y(2)*L-2*lambdana*Lsa+log(R)-delta);
% Change of Population Inversion density
dy(2) = Wp-gamma*sigmaes*c*y(1)*y(2)-(y(2)/tg);
% Change of Photocarrier density
dy(3) = c*y(1)*lambdana-(y(3)/tsa);
end

[darova]

When you put y(1) = 0 as initial condition

% Change of Photon density
dy(1,1) = (y(1)/tr)*(2*sigmaes*y(2)*L-2*lambdana*Lsa+log(R)-delta);

Initial condition y(3) = 0:

% Change of Photocarrier density
dy(3,1) = c*y(1)*lambdana-(y(3)/tsa);

What do you expect to get?

[NUR SIDDIQ]

Thanks for your respond Darova

Should I use random value to generate y(1) ?

[darova]

Can't tell you that. I have no knowledge in this area

%y(1) = Photon Density

What initial Photon Density can be?

[NUR SIDDIQ]

Okay, thank you

It should be start from 1E-10

Answer 1

[Subhadeep Koley]

I suspect the problem is due to the initial conditions. As you stated, if you declare 1e-10 as initial conditions then you will get something in the plot as attached below.

y0=[1e-10; 1e-10; 1e-10];