Lobatto := proc (mu0, alpha, beta, n, x1, x2, x, w)
  #
  # The procedure Lobatto computes the Gauss-Lobatto quadrature rule
  #
  #       b                    n + 2
  #       /                    -----
  #      |                      \
  #      |  omega(x) f(x) dx =   )   w[k] f(x[k])
  #      |                      /
  #     /                      -----
  #     a                      k = 1
  #
  # where the abscissas x1 and x2 on the boundary are prescribed.
  # 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+1
  #     beta[k], k=1..n
  #     n
  #     x1
  #     x2
  # Output:
  #     x[k], k=1..n+2
  #     w[k], k=1..n+2
  #
  #
  # 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).
  #
  # x1 and x2 denote the prescribed abscissas (x1=a and x2=b).
  #
  # x[k] denotes the k. abscissa.
  # w[k] denotes the weight at the abscissa x[k].
  #

  local alphan2tilde, betan1tilde, D, e, i, Q, T, y, z;

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

  # solve (Tn+1 - x1 I) y = en+1 for y[n+1]
  e := linalg[vector] (n+1, 0);  e[n+1] := 1;
  y := linalg[linsolve] (T, e);

  # set up the symmetric tridiagonal matrix T = Tn+1 - x2 I
  for i from 1 to n+1 do
    T[i,i] := alpha[i] - x2;
  od;

  # solve (Tn+1 - x2 I) z = en+1 for z[n+1]
  z := linalg[linsolve] (T, e);

  # compute alphan+2tilde and betan+1tilde
  alphan2tilde := simplify ((x2 * y[n+1] - x1*z[n+1]) /
                            (y[n+1] - z[n+1]));
  betan1tilde := simplify (sqrt ((x2-x1) / 
                            (y[n+1] - z[n+1])));

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

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

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

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

  RETURN (NULL);
end: