Solving laser rate equation with ode45

Hey! I am kinda new to Matlab and I am currently doing my thesis on Semiconductor lasers. I am currently trying to model the rate equations for semiconductor lasers. I have wrote a code but it shows some errors. I have been stuck on this issue for quite some time and it seems that I have done something wrong with ode45. Any type of help is greatly appreciated.

------Main program-----------

ti = 0;
tf = 10;
tspan=[ti tf]; 
y0=[0, 1]; 
% tp = photon lifetime
tp = 1E-12;
% n0 = internal efficiency
n0 = 1E-18;
% a = linewidth broadening factor
a = 2.5E-16;
% L = optical confinement factor
L = 0.3;
% te = electron lifetime
te = 1E-19;
% Rsp = recombination rate
Rsp = ((n0 *1E-18)/te);
% Ng = group refractive index
Ng = 4;
% d = active layer thickness
d = 3E-7;
% g = optical gain
g = 3.16E-12;
b = 2E-5;
J = 2000;
S = 1;

%options? 
[T,Y]= ode45(@rate_eq,tspan,y0,tp,n0,L,Rsp,S,Ng,d,g,b,J); 

plot(T,Y(:,1),'-',T,Y(:,2),'.'); 
title('plot of carrier and photon densities'); 
xlabel('time'); 
ylabel('densities');

-------function--------

function df = rate_eq(tspan,y0,tp,n0,L,Rsp,S,Ng,d,g,b,J,c,e) 

% e = electron charge
e = 1.602E-19;
% c = speed of light
c = 3E8;

% rate equation for carrier density
% Carriers 
df(1,1) = (J/(e * d)) - (Rsp/n0) - ((c*L*g*S)/(Ng));
% Photons 
df(2,1) = ((c*L*g*S)/(Ng)) - (S/tp) + (b*Rsp);

end

Thanks in advance for helping!!

[chandrakanth MS]

can you explain me about the output, like why exactly is it linear? and could give me the refernce from where you got the rate equations from

Answer 1

[Torsten]

ti = 0;
tf = 10;
tspan=[ti tf]; 
y0=[0; 1]; 
% tp = photon lifetime
tp = 1E-12;
% n0 = internal efficiency
n0 = 1E-18;
% a = linewidth broadening factor
a = 2.5E-16;
% L = optical confinement factor
L = 0.3;
% te = electron lifetime
te = 1E-19;
% Rsp = recombination rate
Rsp = ((n0 *1E-18)/te);
% Ng = group refractive index
Ng = 4;
% d = active layer thickness
d = 3E-7;
% g = optical gain
g = 3.16E-12;
b = 2E-5;
J = 2000;
S = 1;
% e = electron charge
e = 1.602E-19;
% c = speed of light
c = 3E8;

[T,Y]= ode45(@(t,y)rate_eq(t,y,tp,n0,L,Rsp,S,Ng,d,g,b,J,c,e),tspan,y0); 
plot(T,Y(:,1),'-',T,Y(:,2),'.'); 
title('plot of carrier and photon densities'); 
xlabel('time'); 
ylabel('densities');

function df = rate_eq(t,y,tp,n0,L,Rsp,S,Ng,d,g,b,J,c,e) 

% rate equation for carrier density
% Carriers 
df(1,1) = (J/(e * d)) - (Rsp/n0) - ((c*L*g*S)/(Ng));
% Photons 
df(2,1) = ((c*L*g*S)/(Ng)) - (S/tp) + (b*Rsp);

end

Best wishes

Torsten.

[Razia Sultana]

ti = 0;
tf = 10;
tspan=[ti tf]; 
y0=[0; 1]; 
% tp = photon lifetime
tp = 1E-12;
% n0 = internal efficiency
n0 = 1E-18;
% a = linewidth broadening factor
a = 2.5E-16;
% L = optical confinement factor
L = 0.3;
% te = electron lifetime
te = 1E-19;
% Rsp = recombination rate
Rsp = ((n0 *1E-18)/te);
% Ng = group refractive index
Ng = 4;
% d = active layer thickness
d = 3E-7;
% g = optical gain
g = 3.16E-12;
b = 2E-5;
J = 2000;
S = 1;
% e = electron charge
e = 1.602E-19;
% c = speed of light
c = 3E8;
[T,Y]= ode45(@(t,y)rate_eq(t,y,tp,n0,L,Rsp,S,Ng,d,g,b,J,c,e),tspan,y0); 
plot(T,Y(:,1),'-',T,Y(:,2),'.'); 
title('plot of carrier and photon densities'); 
xlabel('time'); 
ylabel('densities');
function df = rate_eq(t,y,tp,n0,L,Rsp,S,Ng,d,g,b,J,c,e) 
% this function is not working, showing errors....
%help me out
% rate equation for carrier density
% Carriers 
df(1,1) = (J/(e * d)) - (Rsp/n0) - ((c*L*g*S)/(Ng));
% Photons 
df(2,1) = ((c*L*g*S)/(Ng)) - (S/tp) + (b*Rsp);
end

[Torsten]

I don't get an error message, but solving the equations you implemented is simple:

y1 = ( (J/(e * d)) - (Rsp/n0) - ((c*L*g*S)/(Ng)) ) * t
y2 = 1 + ( ((c*L*g*S)/(Ng)) - (S/tp) + (b*Rsp) )  * t

Answer 2

[Susrutanarayan Chaudhury]

Actually i am unable to solve a coupled laser rate eq in matlab. I am giving these two rate eq and corresponding parameters.

dy/dt=a-y*(x+1+d)...(1)

dx/dt=s+g*x*(y-alpha)....(2)

where x,y the photon intensity and population inversion respectively. The alpha is the modulated loss which is equal toalpha0*(1+m*sin(2*pi*f*t)), where alpha0 is my off set value which I kept 0.2. and 'm' is my modulation index which i kept 0.001. So I have to modulated the loss inside the cavity sinusoidally .

Other parameters are g=1.2E5,

s=5.65;

d=0.3;

a=23.37;

Can any one please help me to solve it in matlab??