%Demo program for polynomial interpolation
clc
clear all
close all
%Data set
n = 6;
x = [0.0; 0.5; 1.0; 6.0; 7.0; 9.0];
y = [0.0; 1.6; 2.0; 2.0; 1.5; 0.0];
%call polyfit to obtain n-1 order fit:
a = polyfit(x,y,n-1)';
disp ('coefficients in decreasing order')
a
%evaluation of interpolant for plotting
xlo=min(x)-0.5;
xhi=max(x)+0.5;
x2 = linspace(xlo,xhi,120)';
f3 = polyval(a,x2);
figure
plot(x,y,'*',x2,f3)
title('Interpolant and data points using polyval')
xlabel ('X')
ylabel ('p(x) and y(x)')
pause
%Evaluation of interpolant for plotting with Horner
z = x2;
m=length(z);
n = length(a);
pval = zeros(m,1);

% coefficients from polyfit are reverse order
% from those in the horner example
a2 = a(n:-1:1);
%to start, set all pval elements = a(n)
pval = a2(n)*(ones(size(z)));
for i = n-1:-1:1
    pval = z.*pval + a2(i);
end
figure
plot(x,y,'*',z,pval)
title('Interpolant and data points using Horner')
xlabel ('x')
ylabel ('p(x) and y(x)')
pause
%enter coefficients for newton divided differences
%a4 = [0; 3.2; -2.4; 0.375757576; -0.052314352; 0.005687776];
ddiff = zeros(n,n);
ddiff(1,:)=y';
for k = 2:n
    for m = 1:n-k+1
        ddiff(k,m)=(ddiff(k-1,m+1)-ddiff(k-1,m))/(x(m+k-1)-x(m));
    end
end
a4(1:n,1)=ddiff(1:n,1);
disp('divided differences to use as coefficients in Newton fit')
a4
f4 = zeros(120,1);
term = ones(120,1);
for k = 1: n
    f4 = f4 + a4(k)*term;
    term = term.*(x2-x(k));
end
figure
plot(x,y,'*',x2,f4)
title('Newton Interpolant and data points')
xlabel ('X')
ylabel ('p(x) and y(x)')
    
