<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def generate_stieltjes( dist, order, accuracy=100, normed=False, retall=False, **kws):
 """ Discretized Stieltjes' method. Args: dist (Dist):
 Distribution defining the space to create weights for. order (int):
 The polynomial order create. accuracy (int):
 The quadrature order of the Clenshaw-Curtis nodes to use at each step, if approximation is used. retall (bool):
 If included, more values are returned Returns: (list):
 List of polynomials, norms of polynomials and three terms coefficients. The list created from the method with ``len(orth) == order+1``. If ``len(dist) > 1``, then each polynomials are multivariate. (numpy.ndarray, numpy.ndarray, numpy.ndarray):
 If ``retall`` is true, also return polynomial norms and the three term coefficients. The norms of the polynomials with ``norms.shape = (dim, order+1)`` where ``dim`` are the number of dimensions in dist. The coefficients have ``shape == (dim, order+1)``. Examples: [q0^2-1.0, q1^2-4.0q1+2.0] [[1. 1. 2.] [1. 1. 4.]] [[0. 0. 0.] [1. 3. 5.]] [[1. 1. 2.] [1. 1. 4.]] q0^2-q0+0.16666667 [[1. 0.0833 0.0056]] """

    from .. import distributions
    assert not distributions.evaluation.get_dependencies(dist)

    if len(dist) > 1:

        # one for each dimension:
        orth, norms, coeff1, coeff2 = zip(*[generate_stieltjes(
            _, order, accuracy, normed, retall=True, **kws) for _ in dist])

        # ensure each polynomial has its own dimension:
        orth = [[chaospy.setdim(_, len(orth)) for _ in poly] for poly in orth]
        orth = [[chaospy.rolldim(_, len(dist)-idx) for _ in poly] for idx, poly in enumerate(orth)]
        orth = [chaospy.poly.base.Poly(_) for _ in zip(*orth)]

        if not retall:
            return orth

        # stack results:
        norms = numpy.vstack(norms)
        coeff1 = numpy.vstack(coeff1)
        coeff2 = numpy.vstack(coeff2)

        return orth, norms, coeff1, coeff2

    try:
        orth, norms, coeff1, coeff2 = _stieltjes_analytical(
            dist, order, normed)
    except NotImplementedError:
        orth, norms, coeff1, coeff2 = _stieltjes_approx(
            dist, order, accuracy, normed, **kws)

    if retall:
        assert not numpy.any(numpy.isnan(coeff1))
        assert not numpy.any(numpy.isnan(coeff2))
        return orth, norms, coeff1, coeff2
    return orth