function [lambda, y, its, norm_r] = potenze (A, x, tol, max_its)
%
% [lambda, y, its, norm_r] = potenze (A, x0, tol, max_its)
%
if ((nargin <= 1) || isempty (x))
  x = randn (size (A, 1), 1);
end
if ((nargin <= 2) || isempty (tol))
  tol = [1e-6, 1e-6];
end
if ((nargin <= 3) || isempty (max_its))
  max_its = 100;
end
if (isscalar (tol))
  tol = [tol, tol];
end
y = x / norm (x);
x = A * y;
its = 1;
lambda = y' * x;
norm_r = norm (x - lambda * y);
while ((norm_r > tol(1) * abs (lambda) + tol(2)) && (its < max_its))
  y =  x / norm (x);
  x = A * y;
  lambda = y' * x;
  its = its + 1;
  norm_r = norm (x - lambda * y);
end
if (norm_r > tol(1) * abs (lambda) + tol(2))
  warning('Impossibile raggiungere la tolleranza richiesta')
end
%!test
%! A = [-2, 1 ,0, 0;1, -2, 1, 0;0, 1, -2 ,1;0, 0, 0, 1];
%! lambda = potenze (A);
%! assert (abs (lambda), max (abs (eig (A))), 1e-5)
