% Perform Cubic Spline interpolation on function func
% evaluated at points x(1:n)
% function [xout,yout,err] = SplineFit(func,x,plt)
%
% Inputs:
%  func = name of function to be interpolated
%  x = array of n points at which to perform interpolation
%  plt = 1 to produce graphics, 0 not o
%
% Outputs:
%  xout = array of 1000 equally spaced points in range of x
%  yout = value of spline interpolant at xout
%  err = maximum  difference |spline(x)-func(x)| at xout
%  if plot = 1, plot of func(xout), yout=spline(xout) 
%               plot of |func(xout)-yout| at 1000 equally space points
%
function [xout,yout,err] = SplineFit(func,x,plt)
xmin = min(x);
xmax = max(x);
xout = xmin:(xmax-xmin)/999:xmax;
n = length(x);
%
%
fxout = eval([func,'(xout)']);
fx = eval([func,'(x)']);
yout = spline(x,fx,xout);
err = max(abs(yout-fxout));
% Plot
if (plt == 1)
   figure(1)
   hold off, clf
   plot(xout,fxout,'b.',xout,yout,'r',x,fx,'bx')
   grid
   title('f(x) = blue, (xi,f(xi)) = blue x, spline(x) = red')
   xlabel(['Maximum error = ',num2str(err)])
   figure(2)
   hold off, clf
   semilogy(xout,abs(fxout-yout),'b'), grid
   title('|f(x) - p(x)|')
   xlabel(['Maximum error = ',num2str(err)])
end