function y=invinterpolacija(f,x0,x1,x2,delta,maxsteps)

% y=INVINTERPOLACIJA(f,x0,x1,x2,delta,maxsteps) dela inverzno kvadratno 
% interpolacijo za funkcijo f z zacetnimi priblizki x0, x1 in x2. Iteracija 
% se konca, ko se zadnja priblizka razlikujeta za manj kot delta ali ko 
% prekoracimo maksimalno stevilo korakov maxsteps. 
%
% Primer uporabe:
% f = @(x) x^3-5*x+1
% invinterpolacija(f,2,2.05,2.1,1e-12,50)

% Bor Plestenjak 2020

if nargin<5, delta=1e-12; end
if nargin<6, maxsteps=50; end

korak = 0;                    						 
while (abs(x2-x1)>delta) && (korak<maxsteps)
    tmp = interpolacija([f(x2) f(x1) f(x0)]',[x2 x1 x0]',0); % vrednost p(0)
    x0 = x1;
    x1 = x2;
    x2 = tmp; 
    korak = korak+1;
    disp(sprintf('%3d:  %s %0.1e',korak,sprintf('%18.15f  ',x2),abs(x2-x1)))  % izpis
    pause(0.1)
end   

y = x2;

% -----------------------------------------------------------------------------------------------

function y=interpolacija(X,Y,x);

% Aitken-Nevillova shema za interpolacijo. Vrne vrednost interpolacijskega 
% polinoma, podanega z točkami X in vrednostmi Y v točki x.

% Bor Plestenjak 2003

n = length(X);
I = zeros(n);
I(:,1) = Y;
for k = 2:n
   for j = k:n
      I(j,k)=1/(X(j)-X(j-k+1))*((x-X(j-k+1))*I(j,k-1)-(x-X(j))*I(j-1,k-1));
   end
end
y = I(n,n);