Approximate the $n$-step transition matrix $P^n$ as $n$ gets very large, given the $4\times 4$ matrix $$P =\pmatrix{0.5 & 0.5 & 0 & 0 \\ 1 & 0 & 0 & 0 & \\ 0 & 0.5 & 0 & 0.5 & \\ 0 & 0 & 0 & 1}$$.
My solution I ran some calculations using Excel, and found that $P^n$ converges to the following matrix: $$\pmatrix{0.666 & 0.333 & 0 & 0 \\ 0.666 & 0.333 & 0 & 0 \\ 0.333 & 0.166 & 0 & 0.5 \\ 0 & 0 & 0 & 1}.$$
My question Is there a systematic way to solve this problem or in general, for any given matrix $P$?
You can use matriz diagonalization; this is representing some matrix $M$ as $M=PDP^{-1}$, where $D$ is a diagonal matrix. The diagonal elements are necessarily the eigenvalues of $M$, and $P$ is a change of basis. In other words, the space on which $M$ acts admits a basis of eigenvalues, and in this basis $M$'s action has a particularly simple representation (simply multiply each eigenvector component by its corresponding eigenvalue).
In this form, it is easy to derive a formula for $M^n$, and you can easily check that we have $M^n=PD^nP^{-1}$ (where calculating $D^n$ is easy precisely because it is diagonal).
At that point, it's just a matter of taking the limit as $n\to \infty$. Wolfram Alpha can calculate an exact expression for the $n$-th power of your matrix (here), and you can check that letting $n\to \infty$ the matrix becomes precisely that which you found empirically, where $0.666$, $0.333$ and $0.166$ are in fact $2/3$, $1/3$ and $1/6$ (respectively).