I have settled this question here; in brief, a FORTRAN implementation 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 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...

I have settled this question here; in brief, a FORTRAN implementation 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 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...

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

I have settled this question here; in brief, a FORTRAN implementation 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 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...