function hermite
nodes=[0.1,0.2,0.3];  % nodes
n=size(nodes);
n=n(2);          % determine the number of nodes 
dub_nodes=zeros(1,2*n); % set up space for dublicated notes
for i = 1:1:n
    dub_nodes(2*i-1)=nodes(i);
    dub_nodes(2*i)=nodes(i);
end
diff=zeros(2*n,2*n); % set up a storage for the divided differences
diff(1,:)=f(dub_nodes); % sets up values of the function at the nodes = zeroth difference
% set up second column of the `` divided differences table'' 
diff(2,2)=f_der(nodes(1));
for i=2:1:n
diff(2,2*i-1)=(diff(1,2*i-1)-diff(1,2*i-2))/(nodes(i)-nodes(i-1));
diff(2,2*i)=f_der(nodes(i));
end
for i=3:1:2*n
for j=i:1:2*n
diff(i,j)=(diff(i-1,j)-diff(i-1,j-1))/(dub_nodes(j)-dub_nodes(j-i+1));
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot the function versus interpolating polynomial 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1);
subplot(1,1,1);
clf; % clear the plot
% Plot the exact function 
x=[-2:.1:1];
plot(x, f(x), '-r');
hold on; % to keep the plot on 
% calculate the approximated function 
m=size(x);
m=m(2);
y=zeros(1,m);
for i=1:1:m
y(i)=Poly(diff,dub_nodes,x(i));
end
plot(x,y,'-b');
hold off;
Poly(diff,dub_nodes,.18)


function P = Poly(diff, nodes, x)
n=size(nodes);
n=n(2);          % determine the number of nodes 
P=0;
for i=1:1:n
    q=1;
    for j=1:1:i-1
    q=q*(x-nodes(j));
    end
    P=P+diff(i,i)*q;
end

function y=f(x);
%y=sin(x);
%y=3.^x
%y=3.*x.*exp(x)-exp(2.*x);
y=x.^2.*cos(x)-3*x;

function y=f_der(x);
%y=cos(x);
%y=ln(3).*3.^x
%y=3.*x.*exp(x)+3.*exp(x)-2.*exp(2.*x);
y=2.*x.*cos(x)-x.^2.*sin(x)-3;