Q:

How to Find Stationary Distributions of a Squared Transition Matrix?

Given a Markov chain with a transition matrix $p$ \begin{pmatrix}
    0 & 0 & .1 & .9 \\
    0 & 0 & .6 & .4 \\
    \
    .8 & .2 & 0 & 0 \\
    .4 & .6 & 0 & 0
\end{pmatrix}
I can easily determine the stationary distribution. However, if I want to find the stationary distributions of $p^2$ \begin{pmatrix}
    .44 & .56 & 0 & 0 \\
    .64 & .36 & 0 & 0 \\
    \
    0 & 0 & .2 & .8 \\
    0 & 0 & .4 & .6
\end{pmatrix}
I have much more difficulty. 
My first step is to set up a system of equations:
$.44π_1+.64π_2+0π_3+0π_4 = π_1$
$.56π_1+.36π_2+0π_3+0π_4 = π_2$
$0π_1+0π_2+.2π_3+.4π_4 = π_3$
$0π_1+0π_2+.8π_3+.6π_4 = π_4$
Since $π_1+π_2+π_3+π_4=1$, I replace any equation in the system with this equation and subtract $π_1, π_2, π_3$ from the other equations:
$-.56π_1+.64π_2+0π_3+0π_4 = 0$
$.56π_1-.64π_2+0π_3+0π_4 = 0$
$0π_1+0π_2-.8π_3+.4π_4 = 0$
$π_1+π_2+π_3+π_4 = 1$
I now have:
\begin{equation}
\begin{pmatrix}
π_1 & π_2 & π_3 & π_4 \\
\end{pmatrix}
*
\begin{pmatrix}
-.56 & .56 & 0 & 1 \\
.64 & -.64 & 0 & 1 \\
    \
0 & 0 & -.8 & 1 \\
0 & 0 & .4 & 1
\end{pmatrix}
=
\begin{pmatrix}
0&0&0&1\\
\end{pmatrix}
\end{equation}
But since this matrix is not invertible, I cannot solve for \begin{pmatrix}
π_1 & π_2 & π_3 & π_4 \\
\end{pmatrix} by multiplying 
\begin{equation}
\begin{pmatrix}
0&0&0&1\\
\end{pmatrix}
*
\begin{pmatrix}
-.56 & .56 & 0 & 1 \\
.64 & -.64 & 0 & 1 \\
    \
0 & 0 & -.8 & 1 \\
0 & 0 & .4 & 1
\end{pmatrix}^{-1}
\end{equation}
Can anyone help me find the stationary distributions for $p^2$?

A:

Your method is rather unusual. There’s no guarantee that the matrix that you form will be invertible. You’ve been lucky so far in only encountering Markov chains that don’t separate into noncommunicating classes, like this one. Each of the classes has its own stationary distribution, and the general solution is a convex combination of them.  
A more conventional way to do this is to directly solve the equation $\mathbf\pi P=\mathbf\pi$, or $\mathbf\pi(P-I)=0$. This can be done by computing the null space of $(P-I)^T$, which in turn can be done via Gaussian elimination. However, since the chain decomposes into two noncommunicating classes, it might be simpler to work out the two individual stationary distributions, call them $\mathbf\pi_A$ and $\mathbf\pi_B$, which you should be able to do using your method on the $2\times2$ submatrices. Any convex combination $(1-\lambda)\mathbf\pi_A+\lambda\mathbf\pi_B$, with $0\le\lambda\le1$, of these two stationary distributions will also be a stationary distribution.  
To continue with your method, proceed as you might for any system of linear equations. First, move the variables to the other of the product by transposing everything, then augment the coefficient matrix with the vector on the right-hand side of the equation, producing $$\pmatrix{-.56&.64&0&0&0\\.56&-.64&0&0&0\\0&0&-.8&.4&0\\1&1&1&1&1}.$$ Perform Gaussian elimination (a.k.a. row.reduction) on this matrix. (Moving the bottom row to the top before doing anything else will make this easier.) There will be an infinite number of solutions, corresponding to the fact mentioned above that there is no unique stationary distribution for this particular chain.

