% Example 9.3
% cubic polymoial fit
% model: y(x) = a0 + a1*x + a2*x^2 + a3*x^3

clear all
close all
clc

x = [1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2];
y = [1.93 1.61 2.27 3.19 3.19 3.71 4.29 4.95 6.07 7.48 8.72 9.34 11.62];
n = length(x);
m = 4; 

A = zeros(m);
b = zeros(m,1);

% Set up the linear system of equations
for i = 1:m
    for j = 1:m
        sumf = 0;
        sumy = 0;
        for k = 1:n
            fvals = f2(x(k));
            sumf = sumf + fvals(i)*fvals(j);
            sumy = sumy + y(k)*fvals(i);
        end
        A(i,j) = sumf;
        b(i) = sumy;
    end
end

% Display the resulting system
disp('A:')
disp(A)
disp('b:')
disp(b)

% Calculate the coefficients
a = A\b;

% Calculate R^2
tmp = a.*f2(x);      
ymodel = sum(tmp);  

sr = sum((y - (ymodel)).^2);
st = sum((y - (1/n).*(sum(y))).^2);
R2 = 1-sr/st;

disp('--------------------------------------------------------------------------')
disp('Model Parameters')
fprintf('      a0        a1         a2         a3         R^2\n')
fprintf('%10.6g %10.6g %10.6g %10.6g %10.4g\n',a(1),a(2),a(3),a(4),R2)


xplot = linspace(min(x),max(x));
tmp = a.*f2(xplot);      % returns a matrix where column i is {a1*f1(xi), a2*f1(xi), a3*f1(xi)...]
yplot = sum(tmp);

figure
hold on
plot(xplot,yplot,'LineWidth',2)
plot(x,y,'ok','MarkerSize',7,'MarkerFaceColor','k')
xlabel('x')
ylabel('y')
xlim([0,2.5])
ylim([0,14])
title('Cubic Polynomial Fit')