x = [5, 10, 15, 20, 25, 30, 35, 40, 45, 50];
y = [17, 24, 31, 33, 37, 37, 40, 40, 42, 41];

% x = [0, 1, 2, 3, 4, 5];
% y = [2.1, 7.7, 13.6, 27.2, 40.9, 61.1];

format long;

plot(x,y,'sr');
xlabel('x');
ylabel('y');

hold all;

[kcroyals, n] = size(x); 

sumx = 0; sumxy = 0; st = 0;
sumy = 0; sumx2 = 0; sr = 0;
order = 2;

for i = 1:n
    sumx = sumx + x(i);
    sumy = sumy + y(i);
    sumxy = sumxy + x(i)*y(i);
    sumx2 = sumx2 + x(i)*x(i);
end

xm = sumx/n;
ym = sumy/n;

for i = 1:order + 1
    for j = 1:i
        k = i + j - 2;
        sum = 0;
        
        for l = 1:n
            sum = sum + (x(l))^k;
        end

        a(i,j) = sum;
        a(j,i) = sum;
    
    end

    sum = 0;
    
    for j = 1:n
        sum = sum + (y(j))*((x(j))^(i-1));
    end
    
    a(i, order+2) = sum;
end

[hei, wid2] = size(a);

wid = wid2 - 1;

baller = zeros(wid, wid);

for o = 1:(wid)
    for p = 1:(wid)
        baller(o,p) = a(o,p);
    end
end

rightballer = zeros(hei, 1);

for p = 1:(hei)
     rightballer(p,1) = a(p, wid+1);    
end

xballer = (baller^(-1))*rightballer;

a0 = xballer(1, 1);
a1 = xballer(2, 1);
a2 = xballer(3, 1);

for i = 1:n
    st = st + (y(i) - ym)^2;
    sr = sr + (y(i) - a0 - a1*x(i) - a2*((x(i))^2))^2;
end

syx = (sr/(n-2))^(0.5);
r2 = (st-sr)/st;

disp('Parabolic fit');
disp(['a0 = ' num2str(a0)]);
disp(['a1 = ' num2str(a1)]);
disp(['a2 = ' num2str(a2)]);
disp(['r^2 = ' num2str(r2)]);

jared = a0 + a1.*x + a2.*(x.^2);

plot(x, jared);


