% SHOWSUPERINTERP  Generates two figures from SUPERINTERP data.
% One figure shows the points used in interpolation.  Other figure
% shows the polynomial doing the extrapolation.

superinterp  % run it;  now "f", "a", "b", "h", "PNPI" are defined

set(0,'defaultlinelinewidth',2.0,'defaultlinemarkersize',14.0)
figure(1)
xp = a:(b-a)/1000:b;
for k = 1:5
  N = 2^(k-1);  x = a:h(k):b;   y = f(x);
  s = ['o' num2str(k) ';N=' num2str(N) ' points;'];  % nested, colored circles
  plot(x,y,s,'markersize',4.0+3*(5-k+1)),  hold on
  s = ['-' num2str(k) ';P_{' num2str(N) '}(x);'];  % curve with same color
  plot(xp,polyval(polyfit(x,y,N),xp),s)
end
yp = f(xp);  plot(xp,yp,'k;f(x)=cos(3x);')
legend('location','south')
plot([pi pi],[-2 1.5],'--k','linewidth',1.0)
text(pi,-2.4,'3.14159','fontsize',9.0,'rotation',90.0)
axis([0 4 -2 1.5]),  xlabel x,  hold off
print -dpdf allpoints.pdf

figure(2),  plot(h,PNPI,'*'),  xlabel h,  ylabel('P_N(3.14159)')
xp = 0:0.01:4;  yp = polyval(polyfit(h,PNPI,4),xp);
hold on,  plot(xp,yp,'k'),  hold off,  grid on
axis([0 4 -2 3]),  print -dpdf extrapolatethis.pdf
