def romberg(function, a, b, args=(), tol=<NUM_LIT>, rtol=<NUM_LIT>, show=False,<EOL>divmax=<NUM_LIT:10>, vec_func=False):
if isinf(a) or isinf(b):<EOL><INDENT>raise ValueError("<STR_LIT>")<EOL><DEDENT>vfunc = vectorize1(function, args, vec_func=vec_func)<EOL>n = <NUM_LIT:1><EOL>interval = [a,b]<EOL>intrange = b-a<EOL>ordsum = _difftrap(vfunc, interval, n)<EOL>result = intrange * ordsum<EOL>resmat = [[result]]<EOL>err = np.inf<EOL>for i in xrange(<NUM_LIT:1>, divmax+<NUM_LIT:1>):<EOL><INDENT>n = n * <NUM_LIT:2><EOL>ordsum = ordsum + _difftrap(vfunc, interval, n)<EOL>resmat.append([])<EOL>resmat[i].append(intrange * ordsum / n)<EOL>for k in range(i):<EOL><INDENT>resmat[i].append(_romberg_diff(resmat[i-<NUM_LIT:1>][k], resmat[i][k], k+<NUM_LIT:1>))<EOL><DEDENT>result = resmat[i][i]<EOL>lastresult = resmat[i-<NUM_LIT:1>][i-<NUM_LIT:1>]<EOL>err = abs(result - lastresult)<EOL>if err < tol or err < rtol*abs(result):<EOL><INDENT>break<EOL><DEDENT><DEDENT>else:<EOL><INDENT>warnings.warn(<EOL>"<STR_LIT>" % (divmax, err),<EOL>AccuracyWarning)<EOL><DEDENT>if show:<EOL><INDENT>_printresmat(vfunc, interval, resmat)<EOL><DEDENT>return result<EOL>
Romberg integration of a callable function or method.

Returns the integral of `function` (a function of one variable)
over the interval (`a`, `b`).

If `show` is 1, the triangular array of the intermediate results
will be printed.  If `vec_func` is True (default is False), then
`function` is assumed to support vector arguments.

Parameters
----------
function : callable
    Function to be integrated.
a : float
    Lower limit of integration.
b : float
    Upper limit of integration.

Returns
-------
results  : float
    Result of the integration.

Other Parameters
----------------
args : tuple, optional
    Extra arguments to pass to function. Each element of `args` will
    be passed as a single argument to `func`. Default is to pass no
    extra arguments.
tol, rtol : float, optional
    The desired absolute and relative tolerances. Defaults are 1.48e-8.
show : bool, optional
    Whether to print the results. Default is False.
divmax : int, optional
    Maximum order of extrapolation. Default is 10.
vec_func : bool, optional
    Whether `func` handles arrays as arguments (i.e whether it is a
    "vector" function). Default is False.

See Also
--------
fixed_quad : Fixed-order Gaussian quadrature.
quad : Adaptive quadrature using QUADPACK.
dblquad : Double integrals.
tplquad : Triple integrals.
romb : Integrators for sampled data.
simps : Integrators for sampled data.
cumtrapz : Cumulative integration for sampled data.
ode : ODE integrator.
odeint : ODE integrator.

References
----------
.. [1] 'Romberg's method' http://en.wikipedia.org/wiki/Romberg%27s_method

Examples
--------
Integrate a gaussian from 0 to 1 and compare to the error function.

>>> from scipy import integrate
>>> from scipy.special import erf
>>> gaussian = lambda x: 1/np.sqrt(np.pi) * np.exp(-x**2)
>>> result = integrate.romberg(gaussian, 0, 1, show=True)
Romberg integration of <function vfunc at ...> from [0, 1]

::

   Steps  StepSize  Results
       1  1.000000  0.385872
       2  0.500000  0.412631  0.421551
       4  0.250000  0.419184  0.421368  0.421356
       8  0.125000  0.420810  0.421352  0.421350  0.421350
      16  0.062500  0.421215  0.421350  0.421350  0.421350  0.421350
      32  0.031250  0.421317  0.421350  0.421350  0.421350  0.421350  0.421350

The final result is 0.421350396475 after 33 function evaluations.

>>> print("%g %g" % (2*result, erf(1)))
0.842701 0.842701