% set up parameters
N = 35 % maximum fortune 
p = 0.5 %probability of winning
q = 1-p
count = 50  % number of games 
fps = 10;   % frames per second for animations 

%set up transition matrix and initial probabilites
P = zeros(N+1,N+1);
for i = 2:N P(i+1,i) = p; end;
for i = 2:N P(i-1,i) = q; end;
P(1,1) = 1;
P(N+1,N+1) = 1;

U0 = ones(N+1,1)/(N+1);
testvect = ones(1,N+1);
checkP = testvect-testvect*P
checkP0 = testvect*U0-1

% show probabilities
V = U0;
U = U0;
for j=1:count;
   U=P*U;
   V=[V,U];
   end;
   
% as a matrixplot   
figure(1);   
surf(V);


% as an animation
figure(2)
U = U0;
for j=1:count;
   U=P*U;
   plot(0:N,U,'.');
   axis([-0.5,N+0.5,0,1]);
   drawnow;
   pause(1/fps);
   end;
   
% show iterated transition matrices
figure(3)
for j=1:count
   mesh(P^j);
   axis tight
   pause(1/fps)
   %pause
end;
   
   
% speed of convergence
figure(4)
eigplot(P);

%sort eigenvalues and eigenvectors
[Q,D] = eig(P);        % Q is the diagonalizer, D is the diagonal matrix 
[L,J] = sort(diag(D)); % J is the permutation of indices to achieve the sort 
Q = Q(:,J);
D = D(J,J);

%create steady state, transient and subdominant projections
Einf = zeros(N+1,N+1); Einf(N,N)=1; Einf(N+1,N+1)=1;
Pinf = Q*Einf/Q;
Ptr = P-Pinf;
Esd1 = zeros(N+1,N+1); Esd1(1,1)=1; 
Esd2 = zeros(N+1,N+1); Esd2(N-1,N-1)=1;
Psd = Q*(L(1)*Esd1+L(N-1)*Esd2)/Q;

%show evolution of actual, steady state, transient and subdominant distributions
figure(5)
U = U0;
Uss = Pinf*U0;
Usd = U0;
for j=1:count
   subplot(2,2,1);
   plot(0:N,U,'.');
   axis([-0.5,N+0.5,0,1]);
   U = P*U;
   title('actual probabilities')
      
   subplot(2,2,2);
   plot(0:N,Uss,'.');
   axis([-0.5,N+0.5,0,1]);
   title('steady state')
      
   subplot(2,2,3);
   plot(0:N,U-Uss,'.');
   axis([-0.5,N+0.5,-0.5,0.5]);
   title('transient')
   
   
   subplot(2,2,4);
   plot(0:N,Usd,'.');
   Usd = Psd*Usd;
   axis([-0.5,N+0.5,-0.5,0.5]);
   title('subdominant')
   pause(1/fps)
   end;