function l = rqrtz(A,M)

%QR演算法求矩陣全部特徵值

%已知矩陣:A

%迭代步數:M

%求得的矩陣特徵值:l

A = hess(A);

for i=1:M

N = size(A);

n = N(1,1);

u = A(n,n);

[q,r]=qr(A-u*eye(n,n));

A = r*q+u*eye(n,n);

l = diag(A);

end

------------------------------------

A=[0 5 0 0 0 0;1 0 4 0 0 0;0 1 0 3 0 0;0 0 1 0 2 0;0 0 0 1 0 1;0 0 0 0 1 0]

A =

0 5 0 0 0 0

1 0 4 0 0 0

0 1 0 3 0 0

0 0 1 0 2 0

0 0 0 1 0 1

0 0 0 0 1 0

>> rqrtz(A,50)

ans =

-3.2030

3.2030

-1.8837

1.8837

-0.6167

0.6167

>> eig(A)

ans =

-3.3243

3.3243

-1.8892

-0.6167

1.8892