function[Q,t,p_lost]=MMnr(lambda, mu,n,r)
% [Q,t]=MM1inf(lambda, mu, n,r)
%lambda - intenzity (rate) of input Poisson flow, inter-arrival time = 1/lambda [s]
%mu - intenzity (rate) of service line, mean service time = 1/mu [s]
%r capacity of buffer - max. queue length
% Q - evolution of the queue
% t - vector of discrete time dt simulation steps t=0:0.1:3600
dt=0.1; T=60*60;
L=0;S=0; %S-served, L-lost
na=0; % active servers
Q=zeros(1,ceil(T/dt));
        for i=1:ceil(T/dt)
            dQ=0;
            if rand < lambda*dt    %arrival of a customer           
                if Q(i)<n+r %free         
                na=min(Q(i)+1,n);dQ=dQ+1; S=S+1; %na increase only for Q<n
                else        %full
                 L=L+1;
                end
            end
            if rand < na*mu*dt && Q(i)>0        %completion of a service
                dQ=dQ-1;na=min(Q(i)-1,n);       %na decrease only for Q<=n
            end
             Q(i+1)=Q(i)+dQ;
        end
        p_lost=L/(L+S);
t=0:dt:T;
