function [L_vect,g_vect,P_vect, peak, energy]=qswitchsim(TQ,r,Lf,numPoints)

a=10^-3;
P0=0.1;

beta=10^-8; %assume small value of spon emission
gth=Lf; %final loss= final gain for threshold[lve
g=r*gth; %initial gain value
P=P0;
Li=Lf*r+Lf*.5;
dL=(Lf-Li)/TQ;
L=Li;
S=1-L;
R=(P+1)*gth;

g_vect(1)=g;
P_vect(1)=P;
L_vect(1)=L;

for t=1:1:numPoints,
    dg=a*(R-(P_vect(t)+1)*g_vect(t));
    dP=((S*(exp(g_vect(t))))-1)*P+beta*g_vect(t);
    g_vect(t+1)=g+dg;
    g=g_vect(t+1);
    P_vect(t+1)=P+dP;
    P=P_vect(t+1);
    if(t<=TQ)
        L_vect(t+1)=L+dL;
        L=L_vect(t+1);
        S=1-L;
    else
        L_vect(t+1)=L;
        S=1-L;
    end
end



subplot(2,1,1)
plot(0:1:numPoints,L_vect,'r')
hold on
plot(0:1:numPoints,g_vect)
title(['r=' num2str(r) ', T_Q=' num2str(TQ), ', L_f=' num2str(Lf)])
ylabel('Gain (blue) & Loss (red)')
subplot(2,1,2)
plot(0:1:numPoints,P_vect)
hold on
title(['peak=' num2str(max(P_vect)) ', energy=' num2str(sum(P_vect))])
ylabel('Power')
xlabel('Normalized Time');
peak=max(P_vect);
energy=sum(P_vect);