functions in romberg.i -

             max_doublings= 20  
 
     is the maximum number of times romberg or simpson will split the  
     integration interval in half.  Default is 20.  

             integral= romberg(function, a, b)  
          or integral= romberg(function, a, b, epsilon)  
 
     returns the integral of FUNCTION(x) from A to B.  If EPSILON is  
     given, Simpson's rule is refined until that fractional accuracy  
     is obtained.  EPSILON defaults to 1.e-6.  
     If the notvector= keyword is supplied and non-zero, then FUNCTION  
     may not be called with a list of x values to return a list of  
     results.  By default, FUNCTION is assumed to be a vector function.  
     If the function is not very smooth, simpson may work better.  

             integral= simpson(function, a, b)  
          or integral= simpson(function, a, b, epsilon)  
 
     returns the integral of FUNCTION(x) from A to B.  If EPSILON is  
     given, Simpson's rule is refined until that fractional accuracy  
     is obtained.  EPSILON defaults to 1.e-6.  
     If the notvector= keyword is supplied and non-zero, then FUNCTION  
     may not be called with a list of x values to return a list of  
     results.  By default, FUNCTION is assumed to be a vector function.  
     If the function is very smooth, romberg may work better.