function ipol()
% Compute the interpolating polynomial of a given function

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

% Number of points
  n = 5;

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

% Compute values of y 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

% Plot function and interpolating polynomial

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

  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

  plot(xv, yv, 'b-', xi, yi, 'r*', xv, pv, 'k-');
  legend('f(x)', 'nodes', 'poly', 'Location', 'Best');

end

% The function to be interpolated

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