function interp_demo(int,fun,fprfun)
%
% This demonstrates the behavior of polynomials that interpolate 
% values of a given function and, optionally, its derivatives on 
% equally spaced interpolation points in a given interval. 
%
% If the arguments are int and fun, then the polynomials interpolate 
% function values; if the arguments are int, fun, and fprfun, then 
% the polynomials interpolate both function and derivative values. 
%
% Inputs: 
% Required:
%  int     = interval [a,b] containing the interpolation points
%  fun     = function to be interpolated (function handle); this 
%	     should accept a vector argument and return the 
%	     corresponding vector of function values. 
% Optional:
%  fprfun  = function for evaluating the derivative (function handle); 
%	     this should accept a vector argument and return the 
%	     corresponding vector of function values. 
% 
% Required M-files: div_dif.m, p_eval.m 
% Optional M-file:  dif_dif_hermite.m (required with fprfun). 
%

a = int(1); 
b = int(2); 
if (a >= b), error('The interval [a,b] must have a<b.'); end 
if ~(nargin == 2|nargin == 3)
  error('The number of arguments must be two or three.'); 
end 
xlow = 2\(b+a)-.8*(b-a);
xhigh = 2\(b+a)+.8*(b-a);

h = 100\(xhigh - xlow);
xx = [xlow:h:xhigh];
fvals = feval(fun,xx); 
if nargin == 3
  fprvals = feval(fprfun,xx); 
end

ylow = min(fvals); 
yhigh = max(fvals); 
temp = 2\(ylow+yhigh) + .6*(yhigh-ylow); 
ylow = 2\(ylow+yhigh) - .6*(yhigh-ylow); 
yhigh = temp; 

figure(1); clf; hold on; 
plot([xlow xhigh],[0 0],'k','LineWidth',1);
plot(xx,fvals,'b','LineWidth',1); 
axis([xlow xhigh ylow yhigh]); 

npts = 1; 
while npts > 0 
  npts = input('Type the number of interpolation points (or -1 to quit):'); 
  if npts <= 0, break; end 
  if nargin == 2
    dgr = npts - 1;
  else
    dgr = 2*npts - 1;
  end 
  if npts > 1
    h = (b-a)/(npts - 1); 
    x = [a:h:b]; 
  else
    x = 2\(a+b); 
  end
  y = feval(fun,x); 
  if nargin == 2 
    y = div_dif(x,y);
    ppval = p_eval(xx,x,y);
  else
    ypr = feval(fprfun,x);
    [d,z] = div_dif_hermite(x,y,ypr);
    ppval = p_eval(xx,z,d);
  end 
  delete(findobj('Type','line','Color','r'));
  delete(findobj('Type','line','MarkerFaceColor','r'));
  figure(1); plot(xx,ppval,'r');
  if nargin == 2
    pval = p_eval(x,x,y);
  else 
    pval = p_eval(x,z,d);
  end 
  plot(x,pval,'o',...
         'MarkerSize',4,'MarkerFaceColor','r',...
         'MarkerEdgeColor','k','LineWidth',1.8);
  title(['Polynomial interpolation: degree ', num2str(dgr)]);
end
