function y = gaussQuad( f, a, b, n )
  %--------------------------------------------------------------------------
  % Jonathan R. Senning <jonathan.senning@gordon.edu>
  % Gordon College
  % February 22, 1999
  %
  % Usage: y = gaussQuad(@f, a, b, n)
  %
  % Returns:
  %	y:	estimate of integral of f(x) from a to b
  %
  % Parameters:
  %	f:	handle of function to integrate
  %	a:	left endpoint
  %	b:	right endpoint
  %	n:	number of nodes (points to evaluate at)
  %
  % Uses gaussian quadrature to estimate the integral of f from a to b using
  % n nodes.  The interval [a,b] is mapped to the interval [-1,1].
  %
  % The change in interval is made by the substitution t = 2*(x-a)/(b-a)-1
  % where x is in [a,b] and t is in [-1,1].  Solving this for x gives
  % x = 0.5*(b-a) * t + 0.5*(b+a).  Thus dx = 0.5*(b-a) * dt.  So
  %
  %      / b                / 1
  %      |            (b-a) |
  %      | f(x) dx  = ----- | f( 0.5*(b-a) * t + 0.5*(b+a) ) dt
  %      |              2   |
  %      / a                /-1
  %
  %--------------------------------------------------------------------------

  % Call grule() to generate the nodes and weights.  Only (n+1)/2 nodes and
  % (n+1)/2 weights are returned, taking advantage of the symmetry of the
  % roots of the Legendre polynomials.
  [r, w] = grule( n );

  % Compute some useful quantities and initialize the sum.  If the number of
  % nodes is odd then the sum is initialized to w((n+1)/2) * f((b+a)/2),
  % accounting for the node r = 0.  If the number of points requested is even
  % then the sum is initialized to zero.
  alpha = 0.5 * ( b - a );
  beta  = 0.5 * ( b + a );

  if ( 2 * floor( n / 2 ) - n == 0 )
    sum = 0.0;
  else
    sum = w((n+1)/2) * f( beta );
  end

  % Compute the sum and return the estimate of the integral
  for i = 1 : floor( n / 2 )
    x0 = -alpha * r(i) + beta;
    x1 =  alpha * r(i) + beta;
    sum = sum + w(i) * ( f( x0 ) + f( x1 ) );
  end

  y = alpha * sum;

end
