function [G,p] = Gmatrix(A,q)
G = zeros(size(A,1));
n = sum(A,2);

for i = 1:size(A,1)
    for j = 1:size(A,2)
        G(i,j) = q/size(A,2) + (A(j,i) * (1 - q))/n(j);
    end
end

[P,D] = eig(G);
PS = sum(P(:,1),1);
p = P(:,1)/PS;