[X]=e722

function [X]=e722
%
%   Find virtual delay distribution of PH/PH/1 queue
%   Compare the results with a diffusion approximation
%
[a,T]=e721a; [b,S]=e721b;
[A]=phaseak(b,S,a,T);
[R]=gph1_sra(A); 
[ma,va]=phasemv(a,T); [ms,vs]=phasemv(b,S);
rho=ms/ma; mu=abs(rho-1);
s2=(1/ma)*(vs+rho*rho*va); md=s2/(2*mu); mdi=1/md; y=[1-rho];
m=50;
for z=1:1:m
   x=1-rho*exp(-mdi*z); y=[y x];
end
%
%  Exact computation
%
[n1,n2]=size(R); e=ones(n1,1); St=-S*e; SB=St*b;
TT=S+R*SB; b1=rho*b; y1=[]; 
for z=0:1:m
   x=1-sum(b1*expm(z*TT)); y1=[y1 x];
end
X=zeros(m+1,3); X(:,1)=(0:1:m)'; X(:,2)=y1'; X(:,3)=y';