n=7;
xd=linspace(0,3,n);yd=cos(1+xd.^2);  % xdata ja ydata
x=linspace(0,3); fun=cos(1+x.^2)     % laskentapisteet ja funktion arvot
c=polyfit(xd,yd,n-1);                % Interpolaatiopolynomin kerroinvektori c
poly=polyval(c,x);                   % c-kerroinvektorin määräämä polynomi
                                     % lasketaan x-pisteissä
plot(x,fun,'b',x,poly,'r',xd,yd,'o');grid

% Käydään Maplessa derivoimassa

M=280000
w=(x-xd(1)).*(x-xd(2)).*(x-xd(2)).*(x-xd(4)).*(x-xd(5)).*(x-xd(6)).*(x-xd(7));
Rarvio=M*w/gamma(n+1);
Rmax=max(Rarvio)
Rmax =

   21.3376
plot(x,Rarvio,x,fun-poly,'r',xd,zeros(size(yd)),'*k');ylim([-2.5 2.5])
shg;grid

Koska 7. derivaatalla on huikea piikki, tulee virhearvioon iso "klappi". Mihinkään järkevään väliin ei voida puristaa (edes väljään), koska 7. derivaatta tulee toisaalta myös nollaksi jossain.