f = inline('1./(1+x.^2)', 'x');

 nvec = [5:25];
  
% least squares 
 err = [];
 for j=1:length(nvec)
    n = nvec(j);
    m = 10*n;

    x = linspace(-5,5,m+1);
    p = polyfit(x,f(x),n);     % warning: polyfit uses Vandermonde matrices...
    x = linspace(-5,5,n+1);
    pn = polyfit(x,f(x),n);
    xx = linspace(-5,5,1000);

    figure(1), clf
    plot(xx,f(xx), 'k-',  'linewidth', 2); hold on
    plot(xx, polyval(p,xx), 'b--', 'linewidth', 2) 
    plot(xx, polyval(pn,xx), 'r--', 'linewidth', 2) 
    fprintf('maximum error:  %10.7e\n', max(abs(f(xx)-polyval(p,xx))));
    err = [err;max(abs(f(xx)-polyval(p,xx))) max(abs(f(xx)-polyval(pn,xx)))];
    axis([-5 5 -.5 2])
    pause
 end
 figure(2), clf
 semilogy(nvec, err(:,2), 'b-', 'linewidth', 2), hold on
 semilogy(nvec, err(:,1), 'r--', 'linewidth', 2)
 legend('p_n = polynomial interpolant', 'p_n = least squares polynomial',2)
 set(gca,'fontsize',16)
 xlabel('n')
 ylabel('max |f(x)-p_n(x)|')