%
%  Polynomial Interpolation
%




x = [ 0; 1; -1; 2; -2];
f = [ -5; -3; -15; 39; -9];

n = size(x(:),1);
z = (-2:0.01:2)';
m = size(z,1);

disp(' Display the polynomial bases used in interpolation ')
disp(' with nodes ')
disp(x')

disp(' Hit return to continue')
pause


% Monomials
M = ones(m,n);
for i = 1:n
  M(:,i) = z.^(i-1);
end
   

plot( z, M(:,1), '-');  hold on
plot( z, M(:,2), ':');  hold on
plot( z, M(:,3), '-.');  hold on
plot( z, M(:,4), '--'); hold on
plot( z, M(:,5), '-'); hold off
xlabel('x')
title('Monomial Basis')
legend('M_1', 'M_2', 'M_3', 'M_4', 'M_5', -1)
%print -deps monomials.eps

disp(' Hit return to continue')
pause

% Lagrange polynomials
L = ones(m,n);
for i = 1:n
    for j = 1:n
        if( i ~= j )
            L(:,i) = L(:,i).*(z-x(j))./(x(i)-x(j));
        end
    end
end
   

plot( z, L(:,1), '-');  hold on
plot( z, L(:,2), ':');  hold on
plot( z, L(:,3), ':');  hold on
plot( z, L(:,4), '--'); hold on
plot( z, L(:,5), '--'); hold off
xlabel('x')
legend('L_1', 'L_2', 'L_3', 'L_4', 'L_5', -1)
title('Lagrange Basis')
legend('M_1', 'M_2', 'M_3', 'M_4', 'M_5', -1)
%print -deps lag-pol.eps

disp(' Hit return to continue')
pause


% Newton polynomials
N = ones(m,n);
for i = 1:n
    for j = 1:i-1
        N(:,i) = N(:,i).*(z-x(j));
    end
end
   

plot( z, N(:,1), '-');  hold on
plot( z, N(:,2), ':');  hold on
plot( z, N(:,3), '-.');  hold on
plot( z, N(:,4), '-'); hold on
plot( z, N(:,5), '--'); hold off
xlabel('x')
legend('N_1', 'N_2', 'N_3', 'N_4', 'N_5', -1)
title('Newton Basis')
%print -deps newt-pol.eps

disp(' Hit return to continue')
pause


plot( z, L*f);  hold on
plot( x, f, '*');  hold off
xlabel('x')
title('Interpolation Polynomial')
%print -deps int-pol.eps