A = [4 2 -1; 0 1 3; 0 0 0.8];
Y = powerit(A,[1 1 1]',100,10^(-4))

A = [-4 2 -1; 0 2 3; 0 0 0.8];
Y = powerit(A,[1 1 1]',100,10^(-4))

A = [10-3i 2 -1; 0 1 3; 0 0 0.8];
Y = powerit(A,[0 1 1]',100,10^(-4))

A = [-4 2 -1; 0 2 3; 0 0 0.8];
Y = inversepowerit(A,[1,1,1]',100,10^(-4))

A = [-4 2 -1; 0 2 3; 0 0 0.8];
Y = inversepowerit(A-2.2*eye(3),[1,1,1]',100,10^(-4))

format long;
A = [-4 2 -1; 0 3.9 3; 0 0 0.8];
Y = rayleigh(A,[-1,-1,1]',100,10^-15);
Y(:,4:5)

Y = rayleigh(A,[1,1,1]',100,10^-15);
Y(:,4:5)

Y = rayleigh(A,[-1,1,0]',100,10^-15);
Y(:,4:5)

format short;
A = [1 2 2 -2; 0 2 0 -3; 3 -6 -4 6; -1 2 2 -3];
eig(A)'

Y = powerit(A,[1,1,1,1]',100,10^(-10));
[n,m] = size(Y);
lambda1 = Y(n,m-1)

v1 = Y(n,1:m-2)'/norm(Y(n,1:m-2));
e1 = [1 0 0 0]';
H = eye(4)-2*((v1-e1)*(v1-e1)')/((v1-e1)'*(v1-e1));
C = H*A*H'

B = C(2:4,2:4);
eig(B)'

W = powerit(B,[1,1,1]',100,10^(-10));
[n,m] = size(W);
lambda2 = W(n,m-1)

u2 = W(n,1:3)';
b = C(1,2:4);
alpha = (b*u2)/(lambda2-lambda1);
v2 = H'*[alpha; u2]

norm(A*v2-lambda2*v2)
