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.