% GOODINTERP  Script for good result from Chebyshev-points interpolation of
% f(x)=1/(1+x^2).
% See L. N. Trefethen, "Spectral Methods in MATLAB", SIAM 2000.
% ELB 5/2/05
% See also: BADINTERP

% choose N and distribute points
N=40;
x=12*cos(pi*(0:N)/N);   % "Chebyshev extreme points"

% fit degree N polynomial to  1/(1+x^2)
f=inline('1./(1+x.^2)','x');
y=f(x);
% following command reports warning if N > 9 because we unnecessarily 
% convert to monomial basis:
p=polyfit(x,y,N);

% plot the result
xx=-12:0.01:12;
yy=polyval(p,xx);
plot(x,y,'o',xx,yy), axis([-12 12 -1 2])