load nfp;
N=length(nfp);
j=(1:N);
X=[j ;nfp];
A=ones(4,N);
A(2,:)=X(1,:);
A(3,:)=X(1,:).^2;
A(4,:)=X(1,:).^3;
Y=X(2,:);
C=inv(A*A');
B=C*A*Y';
plot(X(1,:),X(2,:),'g-*')
hold on
l=B(4)*j.^3+B(3)*j.^2+B(2)*j+B(1);
plot(j,l)
axis([0 N 0 max(nfp)])