function ipol()
% Compute the interpolating polynomial of a given function
% at the nodes of a uniform partition of [a,b]

% The interval [a,b]
  a = -1;
  b =  1;

% Number of points
  n = 13;

% Uniform partition of [a,b]
  xi = linspace(a, b, n);

% Compute values of f at these points
  yi = f(xi);

% Compute coefficients of the interpolating polynomial

  dd(1:n) = yi(1:n);
%
  for i=2:n
    for j=n:-1:i
      dd(j) = ( dd(j) - dd(j-1) ) / ( xi(j) - xi(j+1-i) );
    end
  end

% Now dd(i), i=1,..,n contains the coefficients of the Newton form of the
% interpolationg polynomial
% p(x) = a_0 + a_1 (x - x_0) + ... + a_n(x - x_0) (x - x_1) * ... * (x - x_{n-1})
% in the sense that dd(1) = a_0, dd(2) = a_1, ... dd(n) = a_n

% Estimate maximum interpolation error by comparing f and its
% interpolating polynomial at np points uniformly distributed in [a,b]

  np = 250;
  xv = linspace(a, b, np);
  yv = f(xv);
  pv = zeros(1,np);

% Compute values of the interpolating polynomial

  pv(1:np) = dd(n);

  for i=1:n-1
    pv(1:np) = dd(n-i) + ( xv(1:np) - xi(n-i) ) .* pv(1:np);
  end

% Compute error

  ev = yv-pv;
  maxerr = max(abs(ev));

% Plot function, interpolating polynomial and error

  plot(xv, yv, 'b-', xi, yi, 'r*', xv, pv, 'k-');
  legend('f(x)', 'nodes', 'poly', 'Location', 'Best');
  figure();
  plot(xv, ev, 'b-', xv, 0, 'k-');
  title('Interpolation error', 'FontSize', 14);
  text(0, 2, num2str('max absolute err = %e', maxerr), 'FontSize', 12);
end

% The function to be interpolated

function y = f(x)
% y = 1 ./ (1 + 25*x.^2);
  y = atan(40*x);
% y = exp(x);
end
