% HEERMITE SUBROUTINE. 
% DATA OR FUNCTION INTERPOLATION USING HERMITE POLYNOMIALS

function [HP,Q]=SHerm(p,N,x,y,yp)

% P: Polynomial value at point p.
% Q: Divided Differences Matrix. Diagonal entries are Hermite polynomials
%    coefficients.
% N: # of points to be interpolated.
% x: Abscissa of interpolation points.
% y: Ordinates of interpolation points.
% yp: Derivatives at interpolation points.

% Input data at first two columns
   for I=1:N
      z(2*(I-1) + 1)=x(I);
      z(2*(I-1) + 2)=x(I);
      Q(2*(I-1) + 1,1)=y(I);
      Q(2*(I-1) + 2,1)=y(I);
      Q(2*(I-1) + 2,2)=yp(I);
   end
   
   % Computing rest of entries in second column.
   for I=1:N-1
      Q(2*I+1,2)=(Q(2*(I+1),1)-Q(2*I,1))/(z(2*I+1) -z(2*I));
   end
   
   % New entries for third column and on  (regular divided differences).
   for I=3:2*N
      for j=3:I
         Q(I,j)=(Q(I,j-1) - Q(I-1,j-1))/(z(I)-z(I-j+1));
      end
   end
   format short
      
  % Hermite polynomial construction and evaluation at a point p.
   HP=Q(1,1);
   Prod=1;
   for I=2:2*N
      Prod=Prod.*(p-z(I-1));
      HP=HP+Q(I,I).*Prod;
   end
   %disp([N,p,P]);