[X]=e721

function [X]=e721
%
%   Find queue length distribution of PH/PH/1 queue at 
%   arbitrary time t
%   Compare the results with a diffusion approximation
%
[a,T]=e721a; [b,S]=e721b;
[A]=phaseak(b,S,a,T);
[R]=gph1_sra(A); 
[m,n]=size(R); I=eye(m); C=1/sum(b/(I-R)); 
[ma,va]=phasemv(a,T); [ms,vs]=phasemv(b,S);
rho=ms/ma; y0=1-rho; y=[y0];
c1=C*rho*b; c2=(C/ma)*b; RIS=(R-I)/S; ck=1.0e-6;
RS=R; RT=I; yk=1;
while yk > ck
   yk=sum(c1*RS)+sum(c2*RT*RIS); y=[y yk];
   RS=RS*R; RT=RT*R;
end
%
%  diffusion approximation
%
up=va*((1/ma)^3)+vs*((1/ms)^3); dm=2*((1/ms)-(1/ma));
tau=dm/up; rho_h=exp(-tau); r1=1-rho_h;
[m]=length(y);
ro=rho; y1=[y0];
for i=2:m
  yk=ro*r1; y1=[y1 yk]; ro=ro*rho_h;
end
m=31; x1=y(1,1:m)'; x2=y1(1,1:m)';
X=zeros(m,3); x=ones(1,m); x=cumsum(x)-1; 
X(:,1)=x'; X(:,2)=x1; X(:,3)=x2;