function [xx, ww] = gaussquad(n, a, b)
%GAUSSQUAD  Gaussian quadrature integration.
%
%   [X, W] = GAUSSQUAD(N, A, B) returns the base points X and weight factors
%   W for integrating a function from A to B using an N-point Gaussian
%   quadrature rule which is exact for a polynomial of degree 2*N-1 or
%   lower.
%
%   [X, W] = GAUSSQUAD(N, C) assumes the interval is from 0 to C.
%
%   [X, W] = GAUSSQUAD(N) assumes the interval is from -1 to 1.
%
%   GAUSSQUAD(...) with no output arguments plots the base points X and the
%   weight factors W.

%   This program is based on an anonymous Matlab program found on the
%   MathWorks FTP server ftp://ftp.mathworks.com.  Extended and rewritten
%   for speed.
%
%   According to the original program, concepts for finding the base points
%   and weight factors can be found on page 93 of "Methods of Numerical
%   Integration" by Philip Davis and Philip Rabinowitz.

%   Author:      Peter J. Acklam
%   Time-stamp:  2002-03-03 13:49:10 +0100
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   % Check number of input arguments.
   error(nargchk(1, 3, nargin));

   % Assign default values to missing arguments.
   switch nargin
      case 1                    % GAUSSQUAD(N)
         b = 1;
         a = -1;
      case 2                    % GAUSSQUAD(N, C)
         b = a;
         a = 0;
   end

   u = 1:n-1;
   u = u ./ sqrt(4*u.^2 - 1);

   % Same as A = diag(u, -1) + diag(u, 1), but faster (no adding).
   A = zeros( n, n );
   A( 2 : n+1 : n*(n-1) ) = u;
   A( n+1 : n+1 : n^2-1 ) = u;

   % Find the base points X and weight factors W for the interval [-1,1].
   [v, x] = eig(A);
   [x, k] = sort(diag(x));
   w = 2*v(1,k)'.^2;

   % Linearly transform from [-1,1] to [a,b].
   x = (b-a)/2*x + (a+b)/2;     % Transform base points X.
   w = (b-a)/2*w;               % Adjust weigths.

   % Return output arguments or plot data.
   if nargout
      xx = x;
      ww = w;
   else
      h = plot(x, zeros(1, n), '.', x, w, '.');
      set(h, 'MarkerSize', 12);
      set(gca, 'XLim', [a b]);
      title(sprintf('Gaussian Quadrature with %d points', n));
      xlabel('Base points');
      ylabel('Weight factors');
   end