function test_intepolation()
%
% This function tests the polynomial interpolation 
% programs polyinterp.m and evalpoly.m
%
%

xi = [0 0.2  0.3 0.7 1  2 2.1]';
yi = [1 0.2  -4  -2  -1  0 1.3]';

[a] = polyinterp(xi,yi); % polynomial interpolation

xp=linspace(min(xi),max(xi),1000)';

yp=evalpoly(a,xp);
 
% polynomial interpolation of DATA 
figure(1)
clf
plot(xi,yi,'ro')
hold
plot(xp,yp,'b');
xlabel('x')
legend('data points','polynomial interpolant');

% polynomial interpolation of a function

f=@(x) 1./(1+x.^2); % Runge function

xi = linspace(-3,3,20)';
xp = linspace(-3,3,1000)';
[a] = polyinterp(xi,f(xi)); 
yp=evalpoly(a,xp);

figure(2)
clf
plot(xi,f(xi),'ro')
hold
plot(xp,f(xp),'r')
plot(xp,yp,'b--');
xlabel('x')

% Gauss-Legendre-Lobatto points (unevenly spaced)
[xi]=get_GLL_points_and_D_matrix(19,-3,3)
[a] = polyinterp(xi,f(xi)); 
yp=evalpoly(a,xp);

figure(3)
clf
plot(xi,f(xi),'ro')
hold
plot(xp,f(xp),'r')
plot(xp,yp,'b--');
xlabel('x')



end

