function P4(alpha,n)
% Displays the steady-state distribution associated with n tiles
% and jump parameter alpha.

%Set up Q...
Q = zeros(n,n);
for j=1:n
    if j<n
       Q(j+1,j) = alpha;
       Q(1:j,j) = (1-alpha)/j;
    else
       Q(:,n) = 1/n;
    end
end
[L,U,P] = lu(Q - eye(n,n));

% We know that U is singular, U(n,n) = 0, and U(1:n-1,1:n-1) is
% nonsingular.
% Find z so Uz = 0. If we set z(n) = 1 then we get the rest of z by solving
%          U(1:n-1,1:n-1)*z(1:n-1) = -U(1:n-1,n)
z = zeros(n,1); 
z(n) = 1; 
z(1:n-1) = U(1:n-1,1:n-1)\U(1:n-1,n);

% Normalize (optional)...
f = z*(10000/sum(z));
plot(f)
title(sprintf('The Stationary vector for n = %1d  and   alpha = %10.6f',n,alpha))
shg



    
    