
function stirap;
a=[1 0 0]';  % Initial population distribution
t=[0 20];    % Time domain
height=100;  % Laser intensity
width=1;     % width of laser pulse
delta=0;     % Field  detuning
%
% O.D.E 
% Using MATLAB built in partial derivative solver
[T,Y] = ode45(@hmltn,t,a,[],delta,height,width);
%
% printing pulses
%
omsall=height./(cosh(width.*(T-6)).^2);
ompall=height./(cosh(width.*(T-9)).^2);

figure
subplot(2,1,1);
plot(T,abs(Y).^2);
title('populations');
ylabel('Populations');
legend('|b_1|^2','|b_2|^2','|b_3|^2')

subplot(2,1,2);
plot(T,[ompall omsall]);
title('Field profile');
ylabel('Field intensity')
legend('probe','pump');


return



function b=hmltn(t,a,delta,height,width)
oms=height./(cosh(width.*(t-6)).^2);
omp=height./(cosh(width.*(t-9)).^2);
H=i.*[0   omp   0; ...
      omp delta oms; ...
      0   oms   0];
b=H*a;
return