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

% Code 1: Lagrange Characteristic Polynomials
% A. Meseguer (2007)
 

clear all ; close all 

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

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

 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

% Plotting nodes at y=0 and y =1
 plot(xn,ones(length(xn)),'ok') ; hold on ; grid on
 plot(xn,zeros(length(xn)),'ok','markerfacecolor',[0 0 0]) 
 axis([a b -2 2])
 drawnow

% Press any key to generate successive Lagrange Polynomials
 for ipol = 0:N
  plot(xx,CP(:,ipol+1),'color',rand(3,1)) ; hold on ; drawnow
  text(xn(ipol+1)+0.05,1.1,['L',int2str(ipol+1)])
  ['Press any key for next polynomial']
pause
 end

