  The numerical integration of an analytical function $f(x)$ using a finite set
of equidistant points can be performed by quadrature formulas like the
Newton-Cotes. Unlike Gaussian quadrature formulas however, higher-order
Newton-Cotes formulas are not stable, limiting the usable order of such
formulas. Existing work showed that by the use of orthogonal polynomials,
stable high-order quadrature formulas with equidistant points can be developed.
We improve upon such work by making use of (orthogonal) Gram polynomials and
deriving an iterative algorithm, together allowing us to reduce the
space-complexity of the original algorithm significantly.
