%  Inverting the Vandermonde matrix in order to find the
  %  interpolation polynomial Pn for sin(x) in the interval [1 2]
  %  Pn(x) = c_1 x^n + ..+ c_n x + c_{n+1}.
  %
  clear all
  for ip = 0:3
     N  = 10*3^ip;      % Number of interpolation points
     dx = 1/N;
     x  = [1:dx:2]';                 
  % Find the interpolation polynomial using
  % the Vandermonde matrix
     V  = vander(x);    % Matlab's command to produce the vandermonde matrix
     C  = V\sin(x);     % Warning: use \, not /
  % Plot the results at grid values other than the 
  % ones used for the interpolation.
     y = x + 0.1*dx*rand(size(x));
     subplot(2,2,ip+1);
     plot(y,polyval(C,y)-sin(y));
     xlabel('x');
     ylabel('Pn(x)-f(x)');
     title(['N = ',num2str(N)]);
  end
  shg;