code: def quad_fejer(order, lower=0, upper=1, growth=False, part=None):
      
      order = numpy.asarray(order, dtype=int).flatten()
      lower = numpy.asarray(lower).flatten()
      upper = numpy.asarray(upper).flatten()
  
      dim = max(lower.size, upper.size, order.size)
  
      order = numpy.ones(dim, dtype=int)*order
      lower = numpy.ones(dim)*lower
      upper = numpy.ones(dim)*upper
  
      composite = numpy.array([numpy.arange(2)]*dim)
  
      if growth:
          results = [
              _fejer(numpy.where(order[i] == 0, 0, 2.**(order[i]+1)-2))
              for i in range(dim)
          ]
      else:
          results = [
              _fejer(order[i], composite[i]) for i in range(dim)
          ]
  
      abscis = [_[0] for _ in results]
      weight = [_[1] for _ in results]
  
      abscis = chaospy.quad.combine(abscis, part=part).T
      weight = chaospy.quad.combine(weight, part=part)
  
      abscis = ((upper-lower)*abscis.T + lower).T
      weight = numpy.prod(weight*(upper-lower), -1)
  
      assert len(abscis) == dim
      assert len(weight) == len(abscis.T)
  
      return abscis, weight
docs: Generate the quadrature abscissas and weights in Fejer quadrature.
  
      Example:
          >>> abscissas, weights = quad_fejer(3, 0, 1)
          >>> print(numpy.around(abscissas, 4))
          [[0.0955 0.3455 0.6545 0.9045]]
          >>> print(numpy.around(weights, 4))
          [0.1804 0.2996 0.2996 0.1804]