function [eig,x,k] = my_power(A, tol, x0)

% Input:  A : an n by n matrix;
%        tol: tolerance 
%        x0 : non-zero initial guess
%
% Output: eig: the eigenvalue corresponding to the spectral radius.
%          x : the approximation of a unit eigenvector
%          k: the number of iterations

x = x0/norm(x0);
eig0 = 0; eig=1;k=0;
while ((abs(eig0-eig)) >= tol)
   k = k +1;
   eig0 = eig; 
   y = A*x;
   eig = x'*y;
   x = y/norm(y);
end