Для нахождения стационарного распределения методом Гаусса из данной системы уравнений, можно воспользоваться методом Гаусса-Зейделя.

Программирование этого метода предполагает последовательное итерационное решение системы уравнений до достижения заданной точности. Ниже приведен пример кода на Python, в котором реализован метод Гаусса-Зейделя для данной системы уравнений:

def gauss_seidel(p1, p2, p3, p4, epsilon):
while True:
new_p1 = 0.1967*p1 + 0.4561*p2 + 0.3321*p3 + 0.2982*p4
new_p2 = 0.0750*new_p1 + 0.0553*p2 + 0.1986*p3 + 0.0325*p4
new_p3 = 0.2239*new_p1 + 0.4863*new_p2 + 0.0291*p3 + 0.3382*p4
new_p4 = 0.5044*new_p1 + 0.0023*new_p2 + 0.4402*new_p3 + 0.3311*p4
if abs(new_p1 - p1) < epsilon and abs(new_p2 - p2) < epsilon and abs(new_p3 - p3) < epsilon and abs(new_p4 - p4) < epsilon:
break
p1, p2, p3, p4 = new_p1, new_p2, new_p3, new_p4
return p1, p2, p3, p4
p1, p2, p3, p4 = gauss_seidel(0.25, 0.25, 0.25, 0.25, 0.0001)
print(f"Стационарное распределение: p1 = {p1}, p2 = {p2}, p3 = {p3}, p4 = {p4}")

В данном коде функция gauss_seidel выполняет итерации методом Гаусса-Зейделя до тех пор, пока изменения значений не станут меньше заданной точности epsilon. Начальные значения p1, p2, p3, p4 передаются в функцию, и затем они обновляются на каждой итерации.

После выполнения программы будет выведено стационарное распределение.