K = [1/12 1/30 1/60;1/30 1/30 1/42;1/60 1/42 3/140]
M = [1/30 1/105 1/280;1/105 1/280 1/630;1/280 1/630 1/1260]

[evec,eval]=eig(K,M)

xbar=0:.01:1;
xbar = xbar';
phi1=(1-xbar).^2;
phi2=xbar.*phi1;
phi3=xbar.^2.*phi1;

phi=[phi1 phi2 phi3];

modes=phi*evec;

% plot them mode shapes and the shape functions.  The modes have
% been scaled so they both have about the same magnitude

plot(xbar,modes)
subplot(4,1,1)
plot(xbar,modes)
subplot(4,1,2)
plot(xbar,modes(:,1)/max(abs(modes(:,1)))*max(phi1),xbar,phi1);
subplot(4,1,3)
plot(xbar,modes(:,2)/max(abs(modes(:,2)))*max(phi2),xbar,phi2);
subplot(4,1,4)
plot(xbar,-modes(:,3)/max(abs(modes(:,3)))*max(phi3),xbar,phi3);
