%% Polynomial Interpolation Example

f = @(x) cos(x);          % given function
a = -pi; b = pi;
n = 6;                    % CHANGE THIS AND HIT CTRL-RETURN TO RERUN
x = linspace(a,b,n);      % use n equidistant nodes on [a,b]
y = f(x);                 % values of f at nodes are given

xt = linspace(a,b,1e3);   % for plotting p(x)
d = divdiff(x,y);         % find coefficients of Newton form of p(x)
yt = evnewt(d,x,xt);      % evaluate p at points xt
subplot(2,1,1); plot(x,y,'bo',xt,f(xt),'r',xt,yt,'b'); grid on
legend('given points','f(x)','p(x)')
err = yt - f(xt);
maxerr = max(abs(err));
subplot(2,1,2); plot(x,y*0,'bo',xt,yt-f(xt),'b'); grid on
title(sprintf('max. error = %6.2e',maxerr))