%
%% Comparison of polyfit() and spline()

clearvars
% f = @(x) besselj(1, x)
% f = @sin
f = @(x) 1./(1+25*x.^2); % runge function
xFine = -1.05:1e-3:1.05;
yFine = f(xFine);

xSamples = -1:0.25:1; % create sample points
ySamples = f(xSamples); % calculate y values corresponding to sample points

% Interpolate a polynomial for sample points
fittedPoly = polyfit(xSamples, ySamples, numel(xSamples)-1); % returns
% coefficents of the fitted polynomial
yPoly = polyval(fittedPoly, xFine); 

% Use a polynomial for sample points
ySpline = spline(xSamples, ySamples, xFine);

% Plotting the results
figure()
hold on
plot(xFine, yFine)
plot(xSamples, ySamples, "*")
plot(xFine, yPoly)
plot(xFine, ySpline)
legend('original f(x)', 'sample points', 'polyfitted f(x)', 'spline f(x)')
grid on
grid minor
xlabel('x')
ylabel('y')
title("Comparison of polyfit() and spline()")
hold off

%%
%% Extracting Each Individual Coefficent Array from spline()

VPA_PRECISION = 3;
splineStruct = spline(xSamples, ySamples)
i = 0;
% print polynomial for each sub-interval 
for row = splineStruct.coefs' % foreach-like for loop
    disp("===== " + num2str(i+1) + "th piece =====")
    symPoly = poly2sym(row);
    symPolyPrint = vpa(symPoly, VPA_PRECISION)
    derivSymPoly = vpa(diff(symPoly), VPA_PRECISION) % why not
    intSymPoly = vpa(int(symPoly), VPA_PRECISION) % why not
    i = i + 1;
end