clear all;
x =[20.5  32.7 51.0 73.2  95.7];   % x coordinates of data
y=[765  826  873  942  1032];     % y coordinates of data  

% clear;
% x= 0 + rand(10,1); %Generates a vector of length 10 with values between 0 and 1
% y= -2 + 4*rand(10,1); %Generates a vector of length 10 with values between 
%                       % -2 and -2+4
%                       
% clear;
% x= linspace(0,2,11);
% y= -2*x + 3;  % measured data is exactly on line y = -2x + 3
% y=y+ 0.3*(-1 + 2*rand(size(y)));  % measured data on line y = -2x + 3 but with noise added
                      
m=length(x);
for n=1: m-1  % Find the best fit polynomial of degree n. n<=m-1
    A=zeros(n+1, n+1);
    b=zeros(n+1,1);
    for k=1:n+1
        b(k,1)=sum(y.*x.^(k-1));
    for j=1:n+1
    	A(k,j) = sum(x.^(k+j-2));
    end
    end
    condA = cond(A);

    coeff=A\b;    %The coefficients of the polynomial solve the system A*coeff=b

    xx=linspace(min(x), max(x), 500);  
    yy=zeros(size(xx));
    for k=1:n+1
        yy = yy + coeff(k).*xx.^(k-1);
    end

    figure(n);  plot(x,y,'o',xx,yy);
    grid
    legend('Measured data','Least squares fit');
    title(['Least Squares fit of degree ',num2str(n),'    a0 = ',...
        num2str(coeff(1)),'  a1 = ',num2str(coeff(2))])
    xlabel(['Condition number of A = ',num2str(condA)])
end
