function p = PolynomialInterpolation(f,n)
% p = PolynomialInterpolation(f,n)
% computes and displays a polynomial
% interpolant to n+1 equally spaced data 
% sampled from the function f with a 
% polynomial of degree n.
% The vector p returned is a vector of values
% of the fitting polynomial evaluated on an 
% equally spaced grid of 201 points in [-5,5].
x = linspace(-5,5,201);
y = feval(f,x);
hold on
plot(x,y,'LineWidth',2)
title(sprintf('Polynomial Interpolation: %d data points, degree %d polynomial', n+1, n));
xi = linspace(-5,5,n+1);
%xi = 5*cos(pi*(0:n)/n);         % Chebyshev nodes
fxi = feval(f,xi);
plot(xi,fxi,'xg','LineWidth',2,'MarkerSize',10)
% Here are the two essential commands used for the polynomial interpolant
c = polyfit(xi,fxi,n);       % computes the coefficients
p = polyval(c,x);        % evaluates the polynomial
plot(x,p,'r','LineWidth',2,'MarkerSize',10)
legend('Data function','Data points', 'Polynomial',-1);
hold off