Solution (source code)

= Solution

Condition on $X_m=k$ and apply the <Markov property> and time homogeneity:
$$
p_{ij}(m+n)=\sum_{k\in S}\mathbb P_i(X_m=k)\mathbb P(X_{m+n}=j\mid X_m=k)
=\sum_{k\in S}p_{ik}(m)p_{kj}(n).
$$
Thus the transition matrices satisfy $P(m+n)=P(m)P(n)$; since $P(1)=P$ and $P(0)=I$, induction gives $P(n)=P^n$. Averaging over the initial state gives
$$
\boxed{\mathbb P(X_n=i)=(\mu P^n)_i}.
$$
If $\mu P=\mu$, induction gives $\mu P^n=\mu$, so this probability equals $\mu_i$ for every $n$. This is stationarity of the initial distribution.