function solve_stationary_distribution()


markov_matrix = [ 0.25, 0.75; 0.75, 0.25];

distribution = [1;0];

distribution_next = zeros(length(distribution),1);

tolerance = 10e-08;
while true
    
    for id = 1:length(distribution)
        for idn = 1:length(distribution)
            distribution_next(idn) = distribution_next(idn) + distribution(id)*markov_matrix(id,idn);
        end
    end
    distance = max(abs(distribution - distribution_next));
    fprintf('The current error is %0.10f.\n', distance)
    distribution = distribution_next;
    distribution_next = zeros(length(distribution),1);
    if distance<tolerance
        break
    end
end

end