function [chain,state]=markov(T,n,s0,V);
T = [ 0.85 0.15; 0.05 0.95 ];
s0 = [0.4 0.6];
n=100;
[r c]=size(T);
for k=1:r;
if sum(T(k,:)) ~= 1;
disp('error using markov function')
disp(['row ',num2str(k),' does not sum to one']);
disp(' it sums to :');
disp([ sum(T(k,:)) ]);
disp(['normalizing row ',num2str(k),'']);
T(k,:)=T(k,:)/sum(T(k,:));
end;
end;
V=[1:r];
[v1 v2]=size(V);
if v1 ~= 1 |v2 ~=r
disp('error using markov function');
disp(['state value vector V must be 1 x ',num2str(r),''])
if v2 == 1 &v2 == r;
disp('transposing state valuation vector');
V=V';
else;
return;
end;
end
if s0 < 1 |s0 > r;
disp(['initial state ',num2str(s0),' is out of range']);
disp(['initial state defaulting to 1']);
s0=1;
end;
X=rand(n-1,1);
s=zeros(r,1);
s(s0)=1;
cum=T*triu(ones(size(T)));
for k=1:length(X);
state(:,k)=s;
ppi=[0 s'*cum];
s=((X(k)<=ppi(2:r+1)).*(X(k)>ppi(1:r)))';
end;
chain=V*state;