I have settled this question [here][1]; in brief, a [FORTRAN implementation][2] of algorithms for solving this problem have been published in the Collected Algorithms of the ACM.

----------

At the risk of being redundant, I'd like to mention here some other things I mentioned in the (now expanded) answer of mine at that other site for completeness' sake.

Robert's warning of the Runge phenomenon happening is a good one, and it does happen if your abscissas are perversely distributed (relatedly, the underlying Vandermonde matrix is ill-conditioned); the equispaced case being among the worst-behaved point distributions. Abscissas that are "nicely distributed" (e.g. Legendre, Chebyshev, or any other point distribution which "clusters" near the ends) will generally ensure that you have a quadrature rule that behaves tamely even for large numbers of points. (As an aside, the [`chebfun` project][3] hinges on functions being nicely approximated by interpolating polynomials with abscissas that are transformed Chebyshev polynomial roots/extrema.)

Lastly, whatever you finally settle with, you will want to perform the sanity check of ensuring that all the weights of your quadrature rule are of the same sign; any change of sign in the weights can lead to subtractive cancellation when you employ the quadrature rule, and you wouldn't want that...

  [1]: https://math.stackexchange.com/questions/55687/55742#55742
  [2]: http://www.netlib.org/toms/655
  [3]: http://www2.maths.ox.ac.uk/chebfun/