% plot Lagrange polynomials on [a,b]

a = -1;
b =  1;
N =  7;

% function to interpolate
f = @(x) cos(pi/2*x);
% points for interpolation
Xb = linspace(a,b,N);
Yb  = arrayfun(f,Xb);
% Newton's divided differences
c = newtdd(Xb,Yb,N);


% plot the interpolation points
figure(1); clf;
fs = 24; %font size
plot(Xb,zeros(size(Xb)),'x', ...
     'LineWidth',3,'MarkerSize',10);
grid on;
xlabel('x', 'FontSize', fs); ylabel('y', 'FontSize', fs);
hold on;

% plot the function
plot(Xf,f(Xf),'--',...
     'LineWidth',3,'MarkerSize',10);

pause;

% points for visualization
Xf = linspace(a,b,100);
ColorSet = lines(N);
for i = 1:N,
    % Newton's polynomial based on points 1..i
    Yf = arrayfun(@(x)nest(i-1,c,x,Xb),Xf);
    % plot it
    plot(Xf,Yf,'-',...
         Xb(1:i), Yb(1:i),'o',...
         'LineWidth',3,'Color',ColorSet(i,:));
    if(i<N), pause; end
end
hold off;