% Programming and Numerical Methods
% MSc in Computational and Applied Physics

% Code 3: Hermite first order interpolation
% A. Meseguer (2007)

clear all ; close all 


% As in Code 1, we first set up the Lagrange poly.

 N = 13 ;                   % Number of nodes in = 0,1,...,N
 a = -5 ; b = 5 ;           % Domain: x in [a,b]
 xn = a + [0:N]'*(b-a)/N;   % Spatial grid: [x_0, x_1, ... , x_N]' - N x 1  

 NN = 300 ; xx = linspace(a,b,NN)' ;  % A finer grid for evaluation - 100 x 1

% Lagrange Char. Pol.
 CP = ones(NN,N+1) ;      % CP is full of 1's 
 for ipol = 0:N
   fnorm = 1 ;
   for jprod = 0:N
    if jprod ~= ipol ; fnorm = fnorm*(xn(ipol+1)-xn(jprod+1)) ; end
   end
   
   for jprod = 0:N
    if jprod ~= ipol ; CP(:,ipol+1) = CP(:,ipol+1).*(xx - xn(jprod+1)) ; end
   end
   CP(:,ipol+1) = CP(:,ipol+1)/fnorm ;
 end

% A's and B's polynomials
  A = zeros(NN,N+1) ; B = zeros (NN,N+1) ;

  for ipol = 0:N
    sumlip = 0 ;
    for ksum = 0:N
     if ksum ~= ipol ; sumlip = sumlip + 1/(xn(ipol+1)-xn(ksum+1));end
    end
     A(:,ipol+1) = (1-2*(xx-xn(ipol+1))*sumlip).*(CP(:,ipol+1).^2);
     B(:,ipol+1) = (xx-xn(ipol+1)).*(CP(:,ipol+1).^2);
  end
 
% Now we define the function to be interpolated
  f = exp(-xn.^2); % Exhibits Runge, lets try to diff
  df = -2*xn.*f ;

% First we plot the function in a finer grid
  ff = exp(-xx.^2) ;
  plot(xx,ff,'-b') ; %axis([-5 5 -1.5 1.5]) ; grid on ; hold on
  hold on
  plot(xn,f,'or','markerfacecolor',[1 0 0]) ; grid on ;%

% We plot interpolating polynomial
  HIN = A*f + B*df;
  plot(xx,HIN,'-m')

% Lets have a look at the error
  figure(2) 
  semilogy(xx,abs(HIN-ff),'ok','markerfacecolor',[0 0 0])
  title('ERROR')

