% This m file computes the mean number in the system and mean waiting
% time for the M/PH/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
% if the state probabilities are wanted  for P{n(t) = j} then
% user must specify
%     Jmax  - the maximum value for j 


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

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

pi0 = 1-rho
L = pi0 *a*(I-R)^(-2)*R*e
W = L/lambda

Lq = L - rho
Wq =  W - meanX 
j = 1;
while j < Jmax 
    y(j)  = pi0 * a * R^j * e;
    j = j+1;
end 
y0 = pi0 
y
