Computing an infinite sum of matrix powers seems a bit nasty, but actually reformulating the problem as a transfer matrix (where legendary -> legendary with probability 1) let's you express the state of the system after N loops simply as $A^N x_0 = x_n$. Then, A^N can be readily computed by diagonalizing the matrix so that $A = V diag(\lambda_1, ...) V^T$. The matrix exponential means that different eigenvalues decay at different rates. The largest should be lambda_1 = 1 which will survive in the infinite N limit, and tells you the steady state solution. This machinery is a bit more general than the expectation based equations that the author solves with Gaussian elimination, and is a nice application of SVD/matrix diagonalization.