x = [1,2,4,7];
y = [1.6,2.0,2.5,3.5];

n = 1;
N = 3;

T = zeros(size(x));
T(1) = sum(y)
T(2) = y*x'

S=zeros(size(x,2)-1);
for i=1:size(S);
  S(i) = sum(x.^(i-1));
endfor

A = [S(1),S(2);
     S(2),S(3)]
b = [T(1);T(2)]
a = A\b
polynomial_value = @(arg)( a(2)*arg + a(1));

for i=1:size(x,2)
 polynomial_y(i) = polynomial_value(x(i));
end
 y
 polynomial_y
error = sqrt((sum(y.-polynomial_y))^2)