x1 = [0, 0, 0, 5, 5, 5, 10, 10, 10, 15, 15, 15, 20, 20, 20, 25, 25, 25, 30, 30, 30];
x2 = [0, 10, 20, 0, 10, 20, 0, 10, 20, 0, 10, 20, 0, 10, 20, 0, 10, 20, 0, 10, 20]; 
y = [14.6, 12.9, 11.4, 12.8, 11.3, 10.3, 11.3, 10.1, 8.96, 10.1, 9.03, 8.08, 9.09, 8.17, 7.35, 8.26, 7.46, 6.73, 7.56, 6.85, 6.20];

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

format long;

plot(x1, y,'sr');
plot(x2, y,'sr');
xlabel('T (C)');
ylabel('c = 10 g/L');

hold all;

[kcroyals, n] = size(x); 

order = 2;

for i = 1:order + 1
    for j = 1:i
        k = i + j - 2;
        sum = 0;
        
        for l = 1:n
            sum = sum + (x1(l))*(x2(l));
        end

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

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

for i = 1:n
    st = st + (y(i) - ym)^2;
    sr = sr + (y(i) - a1*x(i) - a0)^2;
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);


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

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

plot(x, jared);

baseketball = a0 + a1.*(8)+a2.*(8.^2);
disp(['The predicted dissolved oxygen concentration for T = 8 C is ' num2str(baseketball)]);


