function  [ y, error ] = seeInterp(func,x,xNodes,varargin)
% SEEINTERP - creates a plot comparing a function to its Lagrange interpolant
%    [ y, error ]  = seeInterp(func,x,xNodes)
%
%    returns the interpolant P(x) and its error f(x)-P(x) 
%
%    [ y, error ] = seeInterp(func,x,xNodes,'showError')
%
%    plots both the solution and the error
%
%    func: function handle to return f=func(x)
%    x:    vector of points to calculate interpolant
%    xNodes: vector of nodes of interpolant
%    y:      P(x) the interpolant
%   
%    this version calls polyinterp for interpolation
   
   N = length(xNodes);  % number of interpolation points
   xNodes = xNodes(:); % force to column vector
   yNodes = func(xNodes);

   y = polyinterp(xNodes,yNodes,x); % do lagrange interpolation
   error = func(x) - y;  % calculate error

   
   figure;
   plot(x,y,'b',xNodes,yNodes,'ro',x,func(x),'r--','LineWidth',[1.5]);
   grid;
   set(gca,'fontweight','bold','fontsize',[14]);
   xlabel('x');
   ylabel('y');
   legend('Interpolant','Nodes','true function','Location','Best');
   title(sprintf('F(x), P(x), N=%d',N));


   if nargin > 3 % show the error plot
      figure;
      plot(x,error,'b','LineWidth',[1.5]);
      grid;
      set(gca,'fontweight','bold','fontsize',[14]);
      xlabel('x');
      ylabel('error');
      title(sprintf('Error: f(x)-P(x), N=%d',N));
   end