function x=mc(x0,P,k)
%mc(x0,P,k) simulates the evolution of a Markov chain.  x0 is either an 
%initial state in {1, ..., n} or an an initial distribution.  P is the nxn 
%transition matrix.  k (optional) specifies the run length: 0 ... k.
if nargin==2
   k=1;
end
crs=cumsum(P,2);                     % Cumulative row sums of P
if length(x0)>1
   y=find(rand(1)<cumsum(x0),1);     % Generate initial state
   x=[y,zeros(1,k)];                 % Preallocate x
else
   x=[x0,zeros(1,k)];
end
u=rand(1,k);                         % uniform random values
for i=1:k
   x(i+1)=find(u(i)<crs(x(i),:),1);  % i to i+1 chain transition
end
% M. Day, August 13, 2014.