Romberg

function [wynik,blad] = ro(f,a,b,epsilon)


R(1,1)=(b-a)/2*(f(a)+f(b))
k=2;
do


h=b-a/(2^(k-1));
j=2;

R(k,1)=1/2*(R(k-1,1)+h);

for i=1 : 2^(k-2)
  suma1=f(a+(2*i-1)*h);
endfor

for j=1 : k
  R(k,j)=R(k,j-1)+((R(k,j-1)-R(k-1,j-1))/(4^(j-1)-1);
endfor


 k++;

until(abs(R(k-1,k-1)-R(k-k)<epsilon)
wynik = R(k,k)
blad =

endfunction