% Plot the Newton basis polynomials q_0(x), ..., q_N(x) 
% for [a,b] = [0,1] with x_j = j/N.

 figure(1), clf
 
 xx = linspace(0,1,300)';          % fine grid for plotting
 N = 5;                            % number of grid points

 x = linspace(0,1,N+1)';           % grid points, x(j+1) = j/N = x_j
 qj = ones(size(xx));              % q_0(x) = 1
 for j=0:N
    plot(xx, qj, 'k-','linewidth',2), hold on
    plot(x,zeros(size(x)), 'k.','markersize',20)
    axis([0 1 -0.1 1.1])
    set(gca,'fontsize',16)
    xlabel('x'), ylabel('q_j(x)'), title('Newton basis polynomials')                                         
    input(' press <return> ')                                                                                
    qj = qj.*(xx-x(j+1));          % q_{j+1}(x) = q_j(x) * (x-x_j)
 end