Gauss := proc (mu0, alpha, beta, n, x, w)
  #
  # The procedure Gauss computes the Gauss quadrature rule
  #
  #       b                      n
  #       /                    -----
  #      |                      \
  #      |  omega(x) f(x) dx =   )   w[k] f(x[k])
  #      |                      /
  #     /                      -----
  #     a                      k = 1
  #
  # The weight function omega(x) and the interval [a,b] need not be
  # specified directly. Rather the calculation is based on the orthogonal
  # polynomials p[k](x) with respect to the weight function omega(x) on
  # the interval [a,b]. These polynomials must be specified by their three
  # term recurrence relationship.
  #
  #
  # Input:
  #     mu0
  #     alpha[k], k=1..n
  #     beta[k], k=1..n-1
  #     n
  # Output:
  #     x[k], k=1..n
  #     w[k], k=1..n
  #
  #
  # mu0 denotes the moment
  # 
  #             b
  #             /
  #            |
  #     mu0 =  |  omega(x) dx
  #            |
  #           /
  #           a
  # 
  # The coefficients
  #     alpha[k], beta[k]
  # define the sequence of orthogonal polynomials p[k](x) by the three
  # term recurrence relationship
  #     x p[k-1] = beta[k-1] p[k-2] + alpha[k] p[k-1] + beta[k] p[k],
  # where beta[0] = 0, p[-1] = 0, p[0] = sqrt(1/mu0).
  #
  # x[k] denotes the k. abscissa.
  # w[k] denotes the weight at the abscissa x[k].
  #

  local D, i, Q, T;

  # set up the symmetric tridiagonal matrix T
  T := array (1..n, 1..n, symmetric, [[0$n]$n]);
  for i from 1 to n do
    T[i,i] := alpha[i];
  od;
  for i from 1 to n-1 do
    T[i,i+1] := beta[i];
    T[i+1,i] := beta[i];
  od;

  # compute the eigenvalue decomposition T = QDQ'
  D := evalf (Eigenvals (T, Q));

  # abscissas
  x := array (1..n, [seq (D[i], i=1..n)]);

  # weights
  w := array (1..n, [seq (evalf (mu0) * Q[1,i]^2, i=1..n)]);

  RETURN (NULL);
end: