docs: def _qnwnorm1(n):
      """
      Compute nodes and weights for quadrature of univariate standard
      normal distribution
  
      Parameters
      ----------
      n : int
          The number of nodes
  
      Returns
      -------
      nodes : np.ndarray(dtype=float)
          An n element array of nodes
  
      nodes : np.ndarray(dtype=float)
          An n element array of weights
  
      Notes
      -----
      Based of original function ``qnwnorm1`` in CompEcon toolbox by
      Miranda and Fackler
  
      References
      ----------
      Miranda, Mario J, and Paul L Fackler. Applied Computational
      Economics and Finance, MIT Press, 2002.
  
      """
      maxit = 100
      pim4 = 1 / np.pi**(0.25)
      m = int(fix((n + 1) / 2))
      nodes = np.zeros(n)
      weights = np.zeros(n)
  
      for i in range(m):
          if i == 0:
              z = np.sqrt(2*n+1) - 1.85575 * ((2 * n + 1)**(-1 / 6.1))
          elif i == 1:
              z = z - 1.14 * (n ** 0.426) / z
          elif i == 2:
              z = 1.86 * z + 0.86 * nodes[0]
          elif i == 3:
              z = 1.91 * z + 0.91 * nodes[1]
          else:
              z = 2 * z + nodes[i-2]
  
          its = 0
  
          while its < maxit:
              its += 1
              p1 = pim4
              p2 = 0
              for j in range(1, n+1):
                  p3 = p2
                  p2 = p1
                  p1 = z * math.sqrt(2.0/j) * p2 - math.sqrt((j - 1.0) / j) * p3
  
              pp = math.sqrt(2 * n) * p2
              z1 = z
              z = z1 - p1/pp
              if abs(z - z1) < 1e-14:
                  break
  
          if its == maxit:
              raise ValueError("Failed to converge in _qnwnorm1")
  
          nodes[n - 1 - i] = z
          nodes[i] = -z
          weights[i] = 2 / (pp*pp)
          weights[n - 1 - i] = weights[i]
  
      weights /= math.sqrt(math.pi)
      nodes = nodes * math.sqrt(2.0)
  
      return nodes, weights