docs: def qnwnorm(n, mu=None, sig2=None, usesqrtm=False):
      """
      Computes nodes and weights for multivariate normal distribution
  
      Parameters
      ----------
      n : int or array_like(float)
          A length-d iterable of the number of nodes in each dimension
  
      mu : scalar or array_like(float), optional(default=zeros(d))
          The means of each dimension of the random variable. If a scalar
          is given, that constant is repeated d times, where d is the
          number of dimensions
  
      sig2 : array_like(float), optional(default=eye(d))
          A d x d array representing the variance-covariance matrix of the
          multivariate normal distribution.
  
      Returns
      -------
      nodes : np.ndarray(dtype=float)
          Quadrature nodes
  
      weights : np.ndarray(dtype=float)
          Weights for quadrature nodes
  
      Notes
      -----
      Based of original function ``qnwnorm`` in CompEcon toolbox by
      Miranda and Fackler
  
      References
      ----------
      Miranda, Mario J, and Paul L Fackler. Applied Computational
      Economics and Finance, MIT Press, 2002.
  
      """
      n = np.atleast_1d(n)
      d = n.size
  
      if mu is None:
          mu = np.zeros(d)
      else:
          mu = np.atleast_1d(mu)
  
      if sig2 is None:
          sig2 = np.eye(d)
      else:
          sig2 = np.atleast_1d(sig2).reshape(d, d)
  
      if all([x.size == 1 for x in [n, mu, sig2]]):
          nodes, weights = _qnwnorm1(n[0])
      else:
          nodes = []
          weights = []
  
          for i in range(d):
              _1d = _qnwnorm1(n[i])
              nodes.append(_1d[0])
              weights.append(_1d[1])
  
          nodes = gridmake(*nodes)
          weights = ckron(*weights[::-1])
  
      if usesqrtm:
          new_sig2 = la.sqrtm(sig2)
      else:  # cholesky
          new_sig2 = la.cholesky(sig2)
  
      if d > 1:
          nodes = nodes.dot(new_sig2) + mu  # Broadcast ok
      else:  # nodes.dot(sig) will not be aligned in scalar case.
          nodes = nodes * new_sig2 + mu
  
      return nodes.squeeze(), weights