xdata = [1 2 3 4 5 6 7 8 9 10];
ydata = [1.71 3.43 5.05 6.81 8.4 10.14 11.99 13.53 15.2 17.05];

%number of data points
mmax = 10;
n=2;

for i=1:n
for j=1:n
M(i,j) =0;
end
b(i) = 0;
end

for i=1:mmax
  M(1,1) = M(1,1) + 1;
  M(1,2) = M(1,2) + xdata(i);
  M(2,1) = M(2,1) + xdata(i);
  M(2,2) = M(2,2) + xdata(i)^2;
  b(1) = b(1) + ydata(i);
  b(2) = b(2) + ydata(i)*xdata(i);
end

tolerance = 0.001;

max_it = 100;
a = [0 0];

 output = [0 a];

a=a';
b=b';

for k=1:max_it

    anew(1) =  (b(1) - M(1,2:n)*a(2:n))/M(1,1);
    anew(n) = (b(n) - M(n,1:n-1)*anew(1:n-1)')/M(n,n);
    norm1 = max(abs(anew-a'));
    norm2 = max(abs(anew));
    if ((norm1/norm2)<tolerance)
        output = [output; k anew];
        break
    end
    output = [output; k anew];
    a=anew';
end

iteration = k

error=0;

for i=1:mmax 
ydata2(i) = a(1) + a(2)*xdata(i);
error = error + (ydata(i) - ydata2(i))^2;
end

error
 
plot(xdata,ydata,'+');
hold on;
plot(xdata,ydata2,'r');
