echo on

N = 200;

p1 = 3*randn([1 N]);
p2 = -0.5*p1+1*randn([1 N]);

clf
plot(p1,p2,'k.')
axis equal
hold on
figure(gcf)
pause

D = [p1; p2];
C = D*D'
[E,L] = eig(C)

diag(C)
diag(L)

sum(diag(C))
sum(diag(L))

% biggest eigenvalue
e1 = E(:,2);
lam1 = L(2,2);

% next biggest
e2 = E(:,1);
lam2 = L(1,1);

s = 5;
plot(s*[0 e1(1)],s*[0 e1(2)],'r','linew',2)
pause
plot(s*[0 e2(1)],s*[0 e2(2)],'b','linew',2)
pause

a1 = e1'*D;
D1 = e1*a1;
plot(D1(1,:),D1(2,:),'cd')
pause

a2 = e2'*D;
D2 = D1 + e2*a2;
plot(D2(1,:),D2(2,:),'go')