% y = romberg(f,x,n,tol)
% Izracun vrednosti integrala z Rombergovo metodo.
% Vhod:
% f : analiticno podana funkcija
% x : integracijski interval [a,b]
% n-1 : maksimalno stevilo delitev intervala
% tol : toleranca
% Izhod:
% y : priblizna vrednost integrala
function y = romberg(f,x,n,tol)
  I = zeros(n);
  h = diff(x);
  a = x;
  y = f(x);
  % trapezna metoda s korakom enakim celemu intervalu
  I(1,1) = h*mean(y);
  for i = 2:n
    % integral izracunamo s polovicnim korakom kot predhodnega
    % vrednosti funkcije racunamo le v novih dodanih vozlih
    % integral izracunamo tako da polovici prejsnjega dodamo vsoto
    % vrednosti funkcije v novih vozlih pomnozeno s korakom h
    h = h/2;
    x = a(1)+h:2*h:a(2);
    s = sum(f(x));
    I(i,1) = I(i-1,1)/2+h*s;
    q = 1;
    % z Rombergovo metodo izracunamo boljse priblizne vrednosti iz prejsnjih
    for j = 2:i
      q = q/4;
      I(i,j) = (I(i,j-1)-q*I(i-1,j-1))/(1-q);
    end;
    % iz zanke izstopimo, ko vrednost I(i,i) zadosca temu pogoju
    if abs(I(i,i)-I(i-1,i-1)) <= tol*abs(I(i,i))
      y = I(i,i);
      return;
    end;
  end;
  y = NaN;
  return;