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))