MATLAB: How to integrate a set of coupled ODE’s

coupled equationlaser rate equationMATLAB and Simulink Student Suiteq-switched laser

Hi everyone
I want to calculate pulse energy (Eout) of a q-switched laser, the equation is below:
Meanwhile phi(t) is depend to 3 coupled equation in laser rate equation [y(1)]. I try to use this code to solve the integral equation:
g = @(t) y(1);
phi = integral (g,0,inf,'ArrayValued',1)
Meanwhile after running the code, Matlab still continue calculating more than 2 hours! Anyone know the problem? What I want to get is single value of phi.
This is my full code, you can run by yourself.
Thank you
clear
clc
ti = 0; 
tf = 12E-6; 
tspan=[ti tf]; 
y0=[0.1; 0; 0];
[t,y] = ode45(@rate_eq,tspan,y0); 
%y(1) = Photon Density
%y(2) = Inverted Population Density
%y(3) = Photocarrier Density
subplot(3,1,1);
plot(t,y(:,1)) 
title('IPhoton Density'); 
xlabel('Time'); 
ylabel('Photon Density');
subplot(3,1,2);
plot(t,y(:,2)); 
title('Inverted Population Density '); 
xlabel('Time'); 
ylabel('Inverted Population Density'); 
subplot(3,1,3)
plot(t,y(:,3)); 
title('Photocarrier Density'); 
xlabel('Time'); 
ylabel('Photocarrier Density'); 
 
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 = 2000E-3;
% 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);
%Pulse Energy
g = @(t) y(1);
phi = integral (g,0,inf,'ArrayValued',1)
end

Best Answer

Your system is ‘stiff’, so use one of the stiff solvers instead of ode45: 
[t,y] = ode15s(@rate_eq,tspan,y0);
With this change, you code ran with: 
Elapsed time is 19.170238 seconds.
(on my machine), and produced this plot:
How to integrate a set of coupled ODE's.png
Related Question