Výpočet newtonova intepolačního polynomu + jeho vypsání

i_newton.m

function v = i_newton(x, y, u)
% ======================================================================
% function v = i_newton(x, y, u)
% Pro hodnotu u (vektor hodnot) urci hodnotu v (vektor hodnot) na zaklade 
% Newtonova tvaru interpolacniho polynomu
% Vstup:
%      x... uzly
%      y... fcni hodnoty v uzlech
%      u... vstupni hodnoty (pokud potřebujete pouze ten polynom, zadejte např 1)
% Vystup:
%      v... interpolovane hodnoty
% ======================================================================
 
% test zadani ==========================================================
n = length(x);
if size(x)~=size(y)
  error('nespravne zadane vstupni hodnoty')
end
 
% vypocet pomernych diferenci =========================================
d = zeros(n);    % matice pomernych diferenci je na pocatku nulova
d(:,1)=y';       % prvni sloupec pomernych diferenci 
for k = 2:n     
  for i = 1:n-k+1
    d(i,k) = (d(i+1,k-1) - d(i,k-1))/(x(i+k-1)-x(i));
  end
end
 
%vypocet hodnoty v Newtonove interpolacnim polynomu ====================
dif = d(1,:);    % zajima nas pouze prvni radek pomernych diferenci
v=dif(n);        % pocatecni hodnota polynomu v bode u
for i=(n-1):-1:1 % sestupny cyklus
    v=v.*(u-x(i))+dif(i);
end
%tisk vysledneho polynomu pozor -- = +
fprintf('P(x) = %s + ', strtrim(rats(d(1,1))));
for i=2:length(d(1,:))
  fprintf('%s*',strtrim(rats(d(1,i))));
  for k = 1:i-1
	fprintf('(x-%s)',strtrim(rats(x(k))));
	end
	if length(d(1,:)) ~= i
	fprintf(' + ');
	end
end