We are given a matrix `P=([.6,.1,.1],[.4,.8,.4],[0,.1,.5])` and we are asked to find its steady state vector. That is we want to find vector x such that Px=x.

Px=x => Px-x=0 or (P-I)x=0 where I is the 3x3 identity matrix `I=([1,0,0],[0,1,0],[0,0,1])` .

`P-I=([-.4,.1,.1],[.4,-.2,.4],[0,.1,-.5])` using basic matrix subtraction. (We subtract corresponding elements.)

So `([-.4,.1,.1],[.4,-.2,.4],[0,.1,-.5])([x_1],[x_2],[x_3])=([0],[0],[0])`

We can solve this by reducing the augmented matrix:

`([-.4,.1,.1,|,0],[.4,-.2,.4,|,0],[0,.1,-.5,|,0])`

In reduced row echelon form we get:

`([1,0,-1.5],[0,1,-5],[0,0,0])`

So we have `x_1-1.5x_3=0,x_2-5x_3=0` and `x_3` is free. So the general solution is `([3/2],[5],[1])` . To eliminate the fraction we can choose `x_3=2` to get basis `w=([3],[10],[2])` . The sum of the elements of w is 15; dividing each entry by 15 yields a probability vector `q=([1"/"5],[2"/"3],[2"/"15])` which is the steady state vector.

`Pq=([.6,.1,.1],[.4,.8,.4],[0,.1,.5])([1"/"5],[2"/"3],[2"/"15])=([1"/"5],[2"/"3],[2"/"15])`

With a computer an easier computational way to find the steady state is to realize that this process gives the percentages—the distribution will tend to the steady state regardless of initial distribution. So if we choose some initial vector, whose entries sum to 1, we can approximate the steady state by taking progressively higher powers of the matrix P.

For example, if we start with d= `([.3],[.3],[.4])` then `P^(10)d~~([.20009766],[.66662822],[.13327412])` which is very close to the actual steady state vector. `P^(25)d` is correct to 9 decimal places.