% CLASS9NOV  Polynomial interpolation explored.

% the data
x = [-1  0  3  5]';
y = [ 1  2 -1  4]';

% the Vandermonde solution, using built-in to generated matrix
van = fliplr(vander(x));
a = van \ y;                                  % get the coeffs for standard form
% plot the result:
figure(1)
xp = -2:0.01:6;                               % a fine grid of x-values at which to plot polynomial
plot(x, y, 'o', ...                           % plot data ...
     xp, a(1)+a(2)*xp+a(3)*xp.^2+a(4)*xp.^3)  %   then polynomial
grid on, axis([-2 6 -8 8])

% the Newton solution, building the matrix "by-hand"
new = [1 0 0 0; 1 1 0 0; 1 4 4*3 0; 1 6 6*5 6*5*2];
c = new \ y;                                  % get the coeffs for Newton's form
% alternative is to use forward substitution code:  c = fs(new,y)
% plot the result:
figure(2)
plot(x, y, 'o', ...
     xp,c(1)+c(2)*(xp-x(1))+c(3)*(xp-x(1)).*(xp-x(2))+c(4)*(xp-x(1)).*(xp-x(2)).*(xp-x(3)))
grid on, axis([-2 6 -8 8])

% Horner's algorithm should be used to evaluate instead
% for example, for Vandermonde version,
%   yp = a(1) + xp .* ( a(2) + xp .* ( a(3) + xp .* a(4) ) );
%   plot(x,y,'o',  xp,yp)

