% Single queue 

mu=4;
lambda=3;
dt=0.01;
nsteps=100000;
Q=zeros(1,nsteps);
Q(1)=0;
for i=2:nsteps
    Q(i)=Q(i-1);
    if rand()< lambda*dt 
        Q(i)=Q(i)+1;
    end;
    if (rand()< mu*dt) && (Q(i)>0) 
        Q(i)=Q(i)-1;
    end;
end;

p=[];
for i=0:max(Q)
    p(end+1)=sum(Q==i)/nsteps;
end;

%%
subplot(2,1,1); cla;
t=(1:nsteps)*dt;
plot(t,Q,'.');
title(sprintf('# people in queue (avg=%g)', mean(Q)));
xlabel('t'); 
ylabel('# people');



%%
subplot(2,1,2); cla;
plot(0:max(Q), p, 'o');
title('Probabilities of having exactly k people in queue');
xlabel('k'); ylabel('\pi_k','rotation',0);
rho=lambda/mu;
k=[0:max(Q)];
hold on;
plot(k, rho.^k*(1-rho),'--r');
legend('simulation', 'exact');
