function y = horner(a,c)

n=length(a)-1;
y=a(n+1);

for k = n:-1:1
    y = a(k)+ c*y;
end

function C = lspoly(x,y,M)

n = length(x);
F = zeros(n,M+1);
for k = 1:M+1
   F(:,k) = x'.^(k-1);
end
A = F'*F;
B = F'*y';
C = AB;

clc
clear all
close all

x = [0, 1, 2, 3, 4, 5];
y = [4, -1, 6, 1, -4, -9];

C = lspoly(x,y,2); % finds the parameters
xx = 0:0.01:5;
yy = zeros(1, length(xx));
for i=1:length(xx)
  yy(i) = horner(C,xx(i));
end

CC = lspoly(x,y,3); % finds the parameters
xxx = 0:0.1:5;
yyy = zeros(1, length(xxx));
for i=1:length(xxx)
  yyy(i) = horner(CC,xxx(i));
end

figure(1)
plot(x, y, 'o',  xx, yy, 'r-', xxx, yyy, 'b-')