function hw7p3

rng(0)
m = 50;
A = triu(randn(m));
A = A + triu(A,1)';
niter = 6;

T = A;
for k = 1:niter
  for i = 1:(m-1)
    for j = (i+1):m
      th = 0.5*atan(2*T(i,j)/(T(j,j)-T(i,i)));
      J  = [cos(th),sin(th);-sin(th),cos(th)];
      T([i,j],:) = J'*T([i,j],:);
      T(:,[i,j]) = T(:,[i,j])*J;
    end
  end
end

subplot(1,2,1)
plot(1:m,sort(eig(A)),'r',...
     1:m,sort(diag(T)),'b')
title('eig(A) and diag(T)')
subplot(1,2,2)
plot(1:m,sort(eig(A)) - sort(diag(T)),'k')
title('eig(A) - diag(T)')

keyboard


return
