function [ANS,ATS,arrivet,departt,ET,NS]=MG1Q(IA,S)
% runs a one server queues using the vector of interarrival times IA  and 
% the vector of service times S  of same length. Begins with 0=the initial number in the system.
% intialization step. clock=clock time. Keeps track of NA and ND=number of arrivals,departures
% by present time. EL=system state=(TA,TD)=time of next arrival, time of next departure

%initialize variables
 clock=0;
 NA=0;               % number of arrivals just after current clock time
 ND=0;              % number of departures just after current clock time
 EL=[IA(1),9999];  % times of next events (arrival, departure)
ET=[];             % times of all events, arrivals and departures 
NS=[];              %NS=number in system (immediately after) times ET. 
%initialize the arrival and departure times of customers.
arrivet=cumsum(IA);   % vector of all arrival times
departt=[];           % vector of all departure times
% only run simulation until last customer arrives.
while (clock<arrivet(length(IA)))
  if (EL(1)<EL(2))
     % next event is an  arrival. advance clock and schedule next arrival
     clock=EL(1);
     NA=NA+1;
     
     if (NA<length(IA))  % if this is not the last arrival, schedule another
     EL(1)=clock+IA(NA+1);
  else 
     EL(1)=9999999;
  end
 
     % if system was empty, schedule a  service time for this customer
    if (NA-ND==1)
       EL(2)=clock+S(NA);
    end

  else
     %next event is a departure. advance clock and schedule a new departure.
     clock=EL(2);
     ND=ND+1;
     % schedule next departure  if number in system>0, otherwise reset EL(2)
     if(NA>ND)
        EL(2)=clock+S(ND+1);
     else EL(2)=99999;
     end
    
      departt=[departt clock];
   
end
  ET=[ET,clock];
  NS=[NS NA-ND];
end
% plots the number in system versus the time.
stairs(ET,NS);
%average time in system
ATS=mean(departt-arrivet(1:length(departt)));
% average number in system
k=length(ET);
ANS=sum(NS(1:k-1).*diff(ET))/ET(k);