function [x, w] = grule( n )
  %-------------------------------------------------------------------------
  % Jonathan R. Senning <jonathan.senning@gordon.edu>
  % Gordon College
  % February 21, 1999
  %
  % Usage: [x, w] = grule(n)
  %
  % Returns
  %	x:	Array of roots.  The roots of the polynomial are symmetric
  %		about the origin so only the non-negative roots are returned.
  %	w:	Array of weights.  w(i) is paired with both x(i) and -x(i).
  %
  % Parameters:
  %	n:	Degree of Legendre polynomial.
  %
  % This function produces the roots and associated weights of the Legendre
  % polynomial of degree n.
  %--------------------------------------------------------------------------

  %--------------------------------------------------------------------------
  % I'm no longer certain but I think this came from "Gaussian Quadrature
  % Formulas" by A.H. Stroud and Don Secrest, Prentice-Hall, 1966.
  % The FORTRAN code presented there was transliterated into C during
  % 1986 and then into an Octave/MATLAB M file in February 1999.
  %--------------------------------------------------------------------------

  m = floor((n+1)/2);
  e1 = n*(n+1);

  for i = 0 : m - 1
    t    = (4*i+3)*pi/(4*n+2);
    x0   = (1.0-(1.0-1.0/n)/(8.0*n*n))*cos(t);
    pkm1 = 1.0;
    pk   = x0;

    for k = 2 : n
      t1 = x0*pk;
      pkp1 = t1-pkm1-(t1-pkm1)/k+t1;
      pkm1 = pk;
      pk = pkp1;
    end

    den  = 1.0-x0*x0;
    d1   = n*(pkm1-x0*pk);
    dpn  = d1/den;
    d2pn = (2.0*x0*dpn-e1*pk)/den;
    d3pn = (4.0*x0*d2pn+(2.0-e1)*dpn)/den;
    d4pn = (6.0*x0*d3pn+(6.0-e1)*d2pn)/den;
    u    = pk/dpn;
    v    = d2pn/dpn;
    h    = -u*(1.0+0.5*u*(v+u*(v*v-d3pn/(3.0*dpn))));
    p    = pk+h*(dpn+0.5*h*(d2pn+h/3.0*(d3pn+0.25*h*d4pn)));
    dp   = dpn+h*(d2pn+0.5*h*(d3pn+h*d4pn/3.0));
    h    = h-p/dp;
    fx   = d1-h*e1*(pk+0.5*h*(dpn+h/3.0*(d2pn+0.25*h*(d3pn+0.2*h*d4pn))));
    x(i+1) = x0+h;
    w(i+1) = 2.0*(1.0-x(i+1)*x(i+1))/(fx*fx);
  end

  if ( m+m > n )
    x(m) = 0.0;
  end

end
