%
%  Find the unique n-th polynomial interpolant
%
%   P(x) = p(1) + p(2) x + p(3) x^2 + ... p(n+1)x^n
%
%  for the function f(x) on the interval (a,b) using
%  (n+1) evenly spaced grid points:
%
%    x_i =  a + i*h ,       h=(b-a)/n
%
%  Found by solving M p = y. (See problem #1 of your homework)
%
%  The vector Y returned is the n-th degree polynomial
%  P evaluated at the vector X using the data (x_i,f(x_i)).
%  
% Input   a    Left endpoint of the interval
%         b    Right endpoint of the interval
%         X    Vector of points where we want to 
%              evaluate the interpolating polynomial
%              (so that we may sample the polynomial
%               between the interpolating nodes)
%         n    n+1 is the total number of equally spaced
%              nodes where we interpolate
%
% Output  Y    The values of the interpolating polynomial
%              at the points specified in vector X.
%
% Local Variables
%         x    the set of nodes created by the value of n
%         y    the function values at the interpolating nodes
%         M    matrix used to calculate the interpolating
%              polynomial
%         p    coefficients of the interpolating polynomial

function Y=Pint(a,b,X,n);

% Create the vector of equally spaced nodes
x=a+[0:n]*(b-a)/n;

% Obtain function values at those nodes
y=f(x)';

% Initialize the matrix for the calculation.
M=[ones(n+1,1)];

% Build the matrix for the calculation.
for i=1:n
  M=[M,(x.^i)'];
end; % end-for-loop

% Solve the system for the coefficients p_0, p_1, ...
% Note that the fliplr command is used so that the
% polynomial evaluation below can be done easily.
p=fliplr((M\y)');

% Evaluate the polynomial at the points in the vector X
Y=polyval(p,X);