% 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