How to program finding a stationary distribution using the Gaussian method?

Question

Mathematics

There is a code:

[p1, p2, p3, p4] = solve(...
'p1 = p1*0.1967 + p2*0.4561 + p3*0.3321 + p4*0.2982',...
'p2 = p1*0.0750 + p2*0.0553 + p3*0.1986 + p4*0.0325',...
'p3 = p1*0.2239 + p2*0.4863 + p3*0.0291 + p4*0.3382',...
'p4 = p1*0.5044 + p2*0.0023 + p3*0.4402 + p4*0.3311',...
'p1 + p2 + p3 + p4 = 1');

How to program the finding of the stationary distribution by the Gaussian method based on these lines?

Answer 1

Rewrite the first 3 lines as

p1*(0.1967-1) + p2*0.4561 + p3*0.3321 + p4*0.2982 = 0
p1*0.0750 + p2*(0.0553-1) + p3*0.1986 + p4*0.0325 = 0
p1*0.2239 + p2*0.4863 + p3*(0.0291-1) + p4*0.3382 = 0

The fourth is not needed - it is their linear combination.
Add a line to them
Get a quadratic system of linear equations. You solve it by the Gauss method in any language you know - and that's it.