
function []=odelas
global tn tr tp r

tn=1e-3;
tr=1e-7
r1=1;
r2=0.9;
tp=tr/(1 - r1*r2);
d=1e-3;
sig=1e-22;
hw=4e-19; %[joules]
ifac = 1/(sig*tp/hw);
nth  = tr/(tp*sig*d) 

r=5 ;tp/tn

tspan = [0 400];  % time span in seconds

  
x0 = [0.5; 0.0];
options=odeset('RelTol',1e-8,'AbsTol',1e-8);
[ttrans,ytrans]=ode45(@eom,[-30 0],[0.05; 1.1],options);


x0=ytrans(end,:)
options=odeset('RelTol',1e-8,'AbsTol',1e-8);
[t,y]=ode45(@eomOn,tspan,x0,options);

figure(1);clf
plot(ttrans,ytrans(:,1),ttrans(end)+t,y(:,1))
plot(t,y(:,1),t,y(:,2))

figure(2);clf
subplot(2,1,1)
plot(t*tp*1000,nth*y(:,2))
xlabel 'time in ms'
ylabel 'population inversion density (1/m^3)'
grid
subplot(2,1,2)
plot(t*tp*1000,ifac*y(:,1))
xlabel 'time in ms'
ylabel 'Intensity (W/m^2)'
grid

function [dydt]=eom(t,y)
global tn tp tr
dydt = zeros(2,1);    % a column vector
dydt=[- y(1) +  (tp/tr) * y(2)* y(1);
      tp/tn * (0 - y(2)) -  y(2) * y(1)];
  
function [dydt]=eomOn(t,y)
global tn tp tr r
dydt = zeros(2,1);
dydt=[- y(1) +  (tp/tr) * y(2)* y(1);
      tp/tn * (r - y(2))  -  y(2) * y(1)];

return


