% to use this, type
%    [x,y,y_in,x_node,y_node] = interp(a,b,n)
% the function takes a left-hand node,a, a right-hand node, b, and
% a number of nodes, n.  It needs a subroutine, f.m  It returns
% x, y = f(x), y_in = polynoial interpolation of f(x),
% x_node = the nodes, and y_node = the function values at the
% nodes.  I make the plots with the sequence of commands:
% >> plot(x,y), hold on, plot(x,y_in,'--')
% >> plot(x_node,y_node,'o'), hold off

function [x,y,y_in,x_node,y_node] = interp(a,b,n)


% the points the polynomial will be defined on:
x = a-1:(b-a)/1000:b+1;
[m1,m2] = size(x);
y_in = zeros(size(x));

for i=1:m2
 y(i) = f(x(i));
end	

x_node = zeros(1,n);
y_node = zeros(1,n);

h = (b-a)/(n-1);
for i=1:n
%uniform nodes:
  x_node(i) = a + (i-1)*h;
% Chebyshev nodes:
  x_node(i) = (a+b)/2 + (b-a)/2*cos((i-1)*pi/(n-1));

  y_node(i) = f(x_node(i));
end
for j=1:n
  prod_j = ones(size(x));
  for k=1:n
    if k ~= j
       prod_j = prod_j.*(x-x_node(k))/(x_node(j)-x_node(k));
    end
  end
  y_in = y_in + y_node(j)*prod_j;
end