% This m file computes the mean number in the system and mean waiting
% time for the PH/PH/1 queue
% The user must specifiy
%     a - the initial probability vector for the PH arrivals
%     T - the generator of the nonabsorbing arrival states
%     b - the initial probability vector for the PH services
%     S - the generator of the nonabsorbing service states
[n,m] = size(T);
et = ones(m,1);
It = eye(m);
T0  = -T*et;
meanY = -a*inv(T)*et 
[l,v] = size(S);
es = ones(v,1);
Is = eye(v);
S0 = -S*es;
meanX = -b*inv(S)*es 
temp1 = T0*a;
temp2 = S0*b;

B0 = T ;
B1 = kron(temp1,b);
B2 = kron(It, S0);
A0 = kron(temp1, Is);
A1  = kron(T,Is) + kron(It,S);
A2 = kron(It, temp2);

A1i = inv(A1);
R1 = -A0 *A1i;
[k,g] = size(R1);
er = ones(g,1);
Ir = eye(g); 

R0 = zeros(g,g);
dr = er'*abs(R1- R0)*er;
Rt = R1;
while dr > 1e-4
       R = -Rt^2*A2*A1i - A0*A1i;
       dr = er'*abs(R- Rt)*er;
       Rt = R;
end;

R  
C = -B2* inv(B0)*B1 + A1 + R*A2; 
D = -B2*inv(B0)*et + inv(Ir-R)*er;    
CT = C';
E = [CT(1:g-1, 1:g)
     D'];
F = [zeros(g-1,1)
        1];
pi1 = (E\F)' 
L = pi1*((inv(Ir-R))^2)*er
W = L*meanY   
