code: def density_orthopoly(x, polybase, order=5, xeval=None):
      # polybase = legendre  #chebyt #hermitenorm#
      # polybase = chebyt
      # polybase = FPoly
      # polybase = ChtPoly
      # polybase = hermite
      # polybase = HPoly
  
      if xeval is None:
          xeval = np.linspace(x.min(), x.max(), 50)
  
      # polys = [legendre(i) for i in range(order)]
      polys = [polybase(i) for i in range(order)]
      # coeffs = [(p(x)*(1-x**2)**(-1/2.)).mean() for p in polys]
      # coeffs = [(p(x)*np.exp(-x*x)).mean() for p in polys]
      coeffs = [(p(x)).mean() for p in polys]
      res = sum(c * p(xeval) for c, p in zip(coeffs, polys))
      # res *= (1-xeval**2)**(-1/2.)
      # res *= np.exp(-xeval**2./2)
      return res, xeval, coeffs, polys
  
imports: from scipy import stats, integrate, special
  import numpy
label_desc: Use image node_id 4 for calling a global function with example usage: density_orthopoly(x, polybase, order, xeval) and returns: res,  xeval,  coeffs,  polys