clc, clear all, close all
D=30
iter=30
for i=1:D,
    for j=1:i,
            A(i,j)=rand;
            A(j,i)=A(i,j);
    end
end

r0=rand(D,1);
b0=sqrt(r0'*r0);
q1=r0/b0;
a1=q1'*A*q1;

r1=A*q1-a1*q1;
b1=sqrt(r1'*r1);
q2=r1/b1;
a2=q2'*A*q2;

r2=A*q2-a2*q2-b1*q1;
b2=sqrt(r2'*r2);
q3=r2./b2;
a3=q3'*A*q3;


qs=zeros(D,iter);
qs(:,1)=q1;
qs(:,2)=q2;
qs(:,3)=q3;

as=zeros(iter,1);
as(1)=a1;
as(2)=a2;
as(3)=a3;

bs=zeros(iter,1);
bs(1)=b0;
bs(2)=b1;
bs(3)=b2;

r2=A*q2-a2*q2-b1*q1;
b2=sqrt(r2'*r2);
q3=r2./b2;
a3=q3'*A*q3;

for i=3:iter,
    r=A*qs(:,i)-as(i)*qs(:,i)-bs(i)*qs(:,i-1);
    bs(i+1)=sqrt(r'*r);
    qs(:,i+1)=r./bs(i+1);
    as(i+1)=qs(:,i+1)'*A*qs(:,i+1);
end
e={};
for i=1:iter,
    qs(:,1:i)'*qs(:,1:i)
    T=qs(:,1:i)'*A*qs(:,1:i);
     for i=1:length(T),
         for j=1:length(T),
            T(j+2:length(T),j)=0;
            T(1:j-2,j)=0;
         end
     end
    a=eig(T);
    e{i}=a
end
hold on
for i=1:iter,
    scatter(ones(length(e{i}),1)*i,e{i})
end
accurate=eig(A)
scatter(ones(size(accurate))*iter+1,accurate,'filled')
hold off
Q=[q1 q2 q3];

EYE=Q'*Q;