function[r]=rom(a,b,n,f)
//RETURNS DEFINITE INTEGRAL OF FUNCTION f ON THE INTERVAL <a,b>
//a:real=lower boundary, b:real=upper boundary,..
//n:integer=number of itterations for trapezoidal rule (higher n = more precise result),..
//f:function=function which must be defined separately 
 for t=0:n-1
   h=(b-a)/2^t
   M(t+1,1)=0.5*h*(f(a)+f(b))
   for i=2:2^t
      M(t+1,1)=M(t+1,1)+h*f(a+(i-1)*h)
   end
 end

 for j=1:n-1
  for i=1:n-1
    if(i>=j) then
      M(i+1,j+1)=((4^j)*M(i+1,j) - M(i,j))/((4^j)-1)
    end
  end
 end

 r=M(n,n);
endfunction