% 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 maximum number of nodes, n.  It needs a subroutine, f.m  It returns
% x, y = f(x)
% y_in = a matrix containing all the polynoial interpolations of f(x),
%             from 2 nodes up to n nodes.
% x_node = a matrix containing all the nodes, and
% y_node = a matrix continaing all the function values at the
% nodes.  I make the plots with the sequence of commands:
%
% >> plot(x,y), hold on, plot(x,y_in(1,:),'--')
% >> plot(x_node(1,:),y_node(1,:),'o'), hold off
%
% This will plot the first interpolant.  Similarly,
%
% >> plot(x,y), hold on, plot(x,y_in(2,:),'--')
% >> plot(x_node(2,:),y_node(2,:),'o'), hold off
%
% will plot the second interpolant.

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

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

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

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

for ii=1:n-1

  h = (b-a)/ii;
  for i=1:ii+1
%uniform nodes:
    x_node(ii,i) = a + (i-1)*h;
    y_node(ii,i) = f(x_node(ii,i));
  end
  for j=1:ii+1
    prod_j = ones(size(x));
    for k=1:ii+1
      if k ~= j
         prod_j = prod_j.*(x-x_node(ii,k))/(x_node(ii,j)-x_node(ii,k));
      end
    end
    y_in(ii,:) = y_in(ii,:) + y_node(ii,j)*prod_j;
  end
end