[X]=e649

function [X]=e649
%
%  Example 6.4.9
%
global lm1 lm2 bt1 bt2 b
r=1/.95; rho=0.95;
lm1=1/(r-.20); lm2=1/0.20;
a=[1 0]; T=[-lm1 lm1; 0 -lm2];
[m1,v1]=phasemv(a,T); 
fprintf('  lamdas: %12.8f  %12.8f \n',lm1,lm2);
b=[.05 .95]; bt1=1/10; bt2=.95/(r-0.5);
fprintf('   betas: %12.8f  %12.8f \n',bt1,bt2);
S=[-bt1 0; 0 -bt2];
[m2,v2]=phasemv(b,S);
fprintf('  m1 = %10.4f  v1 = %10.4f \n',m1,v1);
fprintf('  m2 = %10.4f  v2 = %10.4f \n',m2,v2);
%
%  find the root of generalized Erlang
%
a1=fzero('e649_a',0.5);
fprintf(' root for generlized Erlang = %10.8f \n',a1);
%
%  find the root of hyperexponential
%
a2=fzero('e649_b',0.1);
fprintf(' root for hyperexponential  = %10.8f \n',a2);
aa=[a1 a2]
%  
%  Queue length distribution at arrival epochs
%
X=zeros(101,5); z=ones(1,101); z=cumsum(z)-1; X(:,1)=z';
for i=1:2
   a=aa(i); x=1-a; v=[x]; 
   for j=1:100
   x=x*a; v=[v x];
   end
   X(:,1+i)=v';
end
%
%  Queue length distribution at any time
%
for i=1:2
   a=aa(i); x=1-rho; v=[x]; x=rho*(1-a); v=[v x];
   for j=1:99
   x=x*a; v=[v x];
   end
   X(:,3+i)=v';
end