% SUPERINTERP  Attempt to build a super-interpolator by doing successive
% polynomial interpolation with equal spacing  h, and then extrapolating
% to h=0.  Uses Neville's method (nev.m) for both stages.

xx = pi;
f = @(x) cos(3*x);                   % test it on a smooth fcn
a = 0;  b = 4;

h = zeros(1,5);  PNPI = zeros(1,5);  % will extrapolate this data to h=0
for k = 1:5
  N = 2^(k-1);                       % N = 1,2,4,8,16
  h(k) = (b - a) / N;
  x = a:h(k):b;
  y = f(x);
  PNPI(k) = nev(xx,N,x,y);           % P_N(pi)
end

result = nev(0,4,h,PNPI)             % result = F(0)  if  F(h) = P_N(pi)
err = abs(result - f(pi))
