Kwadratura Romberga - Funkcja-oblicz_Rbg

%dk to dokladnosc inaczej takze epsilon
function [wynik_rbg,blad_teo,numer_przyb] = oblicz_Rbg(f,a,b,dk)
  k=1;
  %Obliczenie R(1,1)
  R(k,1)=((b-a)/2)*(f(a)+f(b));
    h(k)=(b-a);
  %Petla odbywajaca się gdy ||Rn-1,n-1 - Rn,n| < Epsilon|
    do
        k++;
        h(k) = h(k-1)/2;    %Obliczenie Sigmy
    s=0;
        for(i=1:(2^(k-2)))
        s+=f(a+((2*i)-1)*h(k));
        endfor
        R(k,1)=(1/2)*(R(k-1,1)+(h(k-1)*s));
        for(j=2:k)
            R(k,j)=R(k,j-1)+((R(k,j-1)-R(k-1,j-1))/(4^(j-1)-1));
        endfor
        blad_teo =abs(R((k-1),(k-1))-R(k,k));
    until(blad_teo < dk)
    wynik_rbg=R(k,k);
  blad_teo;
  numer_przyb = k;
endfunction