function herm(n) 
%  OSCULATING POLYNOMIAL
%p_eval = interp(x_nodes,divdif_y,x_eval)
%
% This is a function
%         p_eval = interp(x_nodes,divdif_y,x_eval)
% It calculates the Newton divided difference form of 
% the interpolation polynomial of degree m-1, where the 
% nodes are given in x_nodes, m is the length of x_nodes, 
% and the divided differences are given in divdif_y.  
% The 
% points at which the interpolation is to be carried out 
% are given in x_eval; and on exit, p_eval contains the 
% corresponding values of the interpolation polynomial.
%
t = 0:.01:1;

h = 1/n;
x_nodes=0:h:1
x_eval=t;

m = length(x_nodes);
y_values = f(x_nodes)
dy_values = fprime(x_nodes)

%%%
divdif_y = y_values;
m = length(x_nodes);
for i=2:m
    for j=m:-1:i
        divdif_y(j) = (divdif_y(j)-divdif_y(j-1)) ...
                           /(x_nodes(j)-x_nodes(j-i+1));
    end
end
%%%
p_eval = divdif_y(m)*ones(size(x_eval));

for i=m-1:-1:1
    p_eval = divdif_y(i) + (x_eval - x_nodes(i)).*p_eval;
end
%%%
figure(1)
title('Lagrange and Hermite interpolation')
hold all
plot(t,f(t),'-',t,p_eval,'.',x_nodes,y_values,'o',t,hermite(x_nodes,y_values,dy_values,t),'--')
h = legend('sin(4\pi*x)','Lagrange','nodes','Hermite',2);
%set(h,'Interpreter','none')
end

function hermite = hermite(x_nodes,y_values,dy_values,t)
% dy_nodes = values of f'
% t = abcissae where "hermite" is computed
nn=max(size(x_nodes));
mm=max(size(t));
hermite = [];
for j=1:mm
    xx=t(j);
    hxv=0;
    for i=1:nn
        den=1;num=1;
        xn=x_nodes(i);
        derLi=0;
        for k=1:nn
            if (k ~= i)
                num=num*(xx-x_nodes(k));
                arg = xn-x_nodes(k);
                den = den*arg;
                derLi = derLi+1/arg;
            end
        end
        Lix2 = (num/den)^2;
        p = (1-2*(xx-xn)*derLi)*Lix2;
        q = (xx-xn)*Lix2;
        hxv = hxv + (y_values(i)*p+dy_values(i)*q);
    end
    hermite = [hermite,hxv];
end
end


function f= f(x)
    f = sin(4*pi*x);
end
function fprime = fprime(x)
    fprime = 4*pi*cos(4*pi*x);
end

%m = length(x_nodes);
%p_eval = divdif_y(m)*ones(size(x_eval));
%for i=m-1:-1:1
%    p_eval = divdif_y(i) + (x_eval - x_nodes(i)).*p_eval;
%end
