import numpy as np A = np.array( [ [0, 1, 1, 0, 0, 0] , [1, 0, 1, 0, 1, 0] , [1, 1, 0, 1, 0, 0] , [0, 0, 1, 0, 1, 0] , [0, 1, 0, 1, 0, 1] , [0, 0, 0, 0, 1, 0] ]) eigvalues, eigvectors = np.linalg.eigh(A) print(eigvalues) print(eigvectors) l_max_i = eigvalues.argmax() l_max = eigvalues[l_max_i] print(l_max_i, l_max) v_max = eigvectors[l_max_i] def some_norm(M): return np.abs(M).max() B = A for k in range(1, 21): B = B.dot(A) print("some kind of error for k =", k, ":", some_norm(B - l_max**k * np.outer(v_max, v_max)) / l_max**k) print(v_max.argmax())