% This m file computes the mean number in the system and mean waiting
% time for the M/PH/1/K+1  queue
% The user must specifiy
%     lambda  - the mean arrival rate of the Poisson arrivals 
%     a - the initial probability vector for the PH services 
%     T - the generator of the nonabsorbing service states
%     k - the size of the queue (system size is k + 1 )

[n,m] = size(T);
e  = ones(m,1);
I  = eye(m);
meanX = -a*inv(T)*e;

RT = lambda*I - lambda*e*a -T;
R = lambda *inv(RT);

sumR = I; 
for i = 1:1:k
  sumR =  sumR +  R^i;
end

pi0 = 1/(a*(sumR-lambda*R^k*inv(T))*e)
Lsum = 0;
for i = 1:1:k 
 y(i) = (pi0*a*R^i)*e; 
 Lsum = i*y(i) + Lsum;
end
 y(k+1) = (pi0*a*R^k*(-lambda*inv(T)))*e;
L = Lsum + (k+1)*y(k+1) 
W = L/lambda
rhoe  = lambda*meanX*(1 - y(k+1)) 

Lq = L - rhoe 
Wq =  W - meanX 
