function poly_interp_plot
% This function calls poly_interp allowing you to
% experiment with polynomial interpolation of different functions
% using (a) equidistant and (b) Chebyshev interpolation knots
%
% CHOOSE A FUNCTION, INTERVAL, DEGREE OF A POLYNOMIAL:
%f=@(x) sin(x);a=-pi;b=pi; n=10; % smooth, all derivatives bounded
%f=@(x) sqrt(x);a=0;b=2;n=10; % first derivative is infinite at x=0
% Runge function: derivatives grow with the order:
f=@(x) 1./(1+x.^2);a=-5;b=5;n=8; 
% INTERPOLATION
[err1,err2,x1,y1,x2,y2,xx,yy,yy1,yy2]=poly_interp(f,a,b,n);
% RESULTS
hold off
plot(xx,yy,'k',xx,yy1,'b',xx,yy2,'r');
title(['n=',int2str(n),';  err_1=',num2str(err1),'  err_2=',num2str(err2)]);
legend('f(x)','equidist.','Chebyshev');xlabel('x');ylabel('y');axis tight
hold on;
plot(x1,y1,'*b',x2,y2,'*r')
end

function [err1,err2,x1,y1,x2,y2,xx,yy,yy1,yy2]=poly_interp(f,a,b,n)
% polynomial interpolation of function f in the int. [a,b]
% using (n+1) sample points: 
% (a) - equidistant, (b) - Chebyshev's points
%      f should allow vector argument!
% p - the polynom, err1 and err2 - maximal errors in [a,b] 
% for the two cases, estimated using nn equidistant query points:
nn=50*n;xx=linspace(a,b,nn);yy=f(xx);

% equidistant sample points
x1=linspace(a,b,n+1);
y1=f(x1); % f in sample points
p1=polyfit(x1,y1,n); % interpolating polynomial
yy1=polyval(p1,xx);
err1=max(abs(yy1-yy)); % error

% Chebyshev's sample points
t=cos((2*(0:n)+1)*pi/(2*n+2));  x2=0.5*(a+b+(b-a)*t);
y2=f(x2); % f in sample points
p2=polyfit(x2,y2,n); % interpolating polynomial
yy2=polyval(p2,xx);
err2=max(abs(yy2-yy)); % error
end
