function X = MarkovChainGS(n,S,P,X1)
% n = the number of slots to be considered
% S = a row vector containing possible states (usually 1:N)
% P = transition probability matrix
% X1 = initial state for slot 1

N = length(S);      % Number of possible states
T = zeros(1,n);     % Preallocation
T(1) = find(S==X1); % Express the states using indices from 1 to N
                    % instead of the provided support S
for k = 2:n
    T(k) =  randsrc(1,1,[S;P(T(k-1),:)]);
end
X = S(T);           % Express the states using the provided support
end