double integrate( string fname, double a, double b, int _n,
  matrix p1, matrix p2, matrix p3, matrix p4, matrix p5)
{
/*
 * Trapezoid rule integration of function named by fname,
 * using _n points over the interval [a,b] with one step of
 * Richardson extrapolation.
 *
 * The function to be integrated must return a row vector
 * of function values when called with a row vector argument.
 *
 * If _n is not specified, 100 points are used by default.
 *
 * p1..p5 are up to 5 extra parameters to be passed directly
 * to the function being integrated.
 *
 * R. Perry, June 1998
 */ 

  int n;
  double s;
  double ans1, ans2;
  double dx;

  if( undef(_n))
    n = 100;
  else
    n = _n;

/*
 * initialize sum using endpoints
 */
  s = sum( feval( fname, [a b], p1, p2, p3, p4, p5))/2;

/*
 * trapezoid rule integration
 */
  dx = (b-a)/n;

  s += sum( feval( fname, a + [1:n-1]*dx, p1, p2, p3, p4, p5));

  ans1 = s*dx;

/*
 * add points between those already computed
 */
  a += dx/2;

  s += sum( feval( fname, a + [0:n-1]*dx, p1, p2, p3, p4, p5));

  ans2 = s*dx/2;

  return (4*ans2-ans1)/3;	/* Richardson extrapolation */
}

/*

Algorithm:

To numerically integrate a function on the interval [a,b], evaluate the
function at n+1 equally spaced points x0, x1, ..., xn (where x0 = a, xn = b,
dx = (b-a)/n) to get the function values f0, f1, ..., fn.

Approximate the area as:

		f0 + f1		f1 + f2			   fn-1 + fn
	area =  ------- dx   +  ------- dx  +	...	+  --------- dx
		   2		   2			        2

The formula can be simplified:

	area =  dx * (  (f0+fn)/2   +  f1  +  f2  +  ...  +  fn-1 )

and implemented in C as:

	s = (f(a) + f(b)) / 2;

	dx = (b-a)/n;

	for( j = 1; j < n; ++j)
	    s += f( a + j*dx );

	ans1 = s * dx;

To get a better estimate of the area, add in more points between those already
computed:
 
	a += dx/2;

	for( j = 0; j < n; ++j)
	    s += f( a + j*dx );

	ans2 = s * dx / 2;

In ans1, dx is the distance between the evaluated points.  In ans2, dx/2
is the distance between the points.  The error in trapezoid rule integration
is mainly proportional to the square of the distance between the points, so:

	ans1 = area + c * dx * dx

	ans2 = area + c * (dx/2) * (dx/2)

If we take  4 * ans2  -  ans1,  the error terms cancel and we are left
with just 3*area.  This is called Richardson extrapolation.

So the final estimate of the area, more accurate than either ans1 or ans2, is:

	(4 * ans2 - ans1) / 3
*/
