The procedure

The Romberg Quadrature is a simple but powerful method to compute numerically the definite integral of functions. The following implementation is taken from the book of Wilhelm Werner (ref. [3]). We made two minor changes:
(1) As initial value we compute instead of R[0,0]:=evalf((limit(f(x),x=a,right)+limit(f(x),x=b,left))*h/2) the value R[0,0]:=evalf((f(a)+f(b))*h/2) which is better suited if the integrand has a pole at one of the endpoints (set f(a)=0 or f(b)=0 in this case).
(2) We made the procedure return the final approximation.