f = @(x) 1./(1+ x.^2); %Function
n= 15; %Number of nodes
x = linspace(-5,5,n); %Nodes
y = f(x); % y data

% Divided differencing algorithm
a = y;
for j=2:n
    for i=n:-1:j
        a(i) = (a(i) - a(i-1))/(x(i) - x(i-j+1));
    end
end

% Create the polynomial and evaluate it
xx = -10:.01:10;
Q = zeros(1,length(x));

for i=1:length(xx)
    Q(i)= evalpol(a,x,n,xx(i));
end

% plot
plot(xx,f(xx),xx,Q,x,y,'o')
axis([-5 5 -1 2])

% This is the function that constructs the polynomial from the divided
% different coefficients a and the nodes x
function Q = evalpol(a,x,n,xx)
    Q = a(n);
    for i=(n-1):-1:1
        Q = Q.*(xx-x(i)) + a(i);
    end
end



