a = -1;
b =  1;

function y = rung( x )
    y = 1 ./ ( 1 + 25 * x .^ 2 );
endfunction

%% equispaced points
n = input( 'polynomial degree = ' );

for i = 1 : n - 1
    h  = ( b - a ) / i;
    xi = a : h : b;   % or   xi = linspace( -1, 1, i );
    yi = rung( xi );

    % polynomial interpolator on equispaced nodes
    x     = linspace( a, b, 321 );
    pequi = polyfit( xi, yi, i );
    g     = rung( x );
    yeq   = polyval( pequi, x );

    %plot( xi, yi, 'o', x, yeq, 'b-', x, g, 'r-' );
    %plot(              x ,yeq, 'b-', x, g, 'r-' );

    err3 = ( g - yeq );
    err4 = norm( err3, 2 );

    plot( x, yeq, 'b-', x, g, 'm-', x, err3, 'r-' );
endfor