function code_plot_Figure4
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Matlab code to generate Figure 4
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


figure
hold on;
grid on;
do_plot(1,1,1)
do_plot(5,1,1) 
do_plot(10,1,1) 

end

function do_plot(b,l,mu)


ind=1;
v_a=0.01:0.01:1000;
for a=v_a
   vpi=compute_stat_dist(a,b,l,mu);
   A=(vpi(4)+(a*vpi(2))/(l+mu+a))/(l+mu+b-(a*b)/(l+mu+a));
   N=[l+a, -mu,-b, 0;
   -l ,mu+a, 0,-b;
   -a ,  0,l+b, 0;
   0,   -a, -l,mu+b];
   v=[vpi(1)+mu*A,vpi(2),vpi(3),vpi(4)]';
   res(ind)=sum(N\v);ind=ind+1;
end
plot(v_a,res)
end





function pi_=compute_stat_dist(a,b,l,mu)

M=[a+l,-mu , -b,-mu; 
   -l ,mu+a,  0,-b; 
   -a ,   0,l+b, 0;
   1  ,   1,  1, 1];
b=[0,0,0,1]'; pi_=M\b;

%M=[-a - l, b, 0, mu, mu, 0 ; ...
%    a, -a - b - l, b, 0, 0, mu ; ...
%    0, a, -b - l, 0, 0, 0 ; ...
%    l, 0, 0, -a - mu, b, 0  ; ...
%    0, l, 0, a, -a - b - mu, b; ...
%    0, 0, l, 0, a, -b - mu ]./(a+b+l+mu)+eye(6);
%mat=M^(100);
%pi_=mat(:,1);
end