docs: def _qnwsimp1(n, a, b):
      """
      Compute univariate Simpson quadrature nodes and weights
  
      Parameters
      ----------
      n : int
          The number of nodes
  
      a : int
          The lower endpoint
  
      b : int
          The upper endpoint
  
      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 ``qnwsimp1`` in CompEcon toolbox by
      Miranda and Fackler
  
      References
      ----------
      Miranda, Mario J, and Paul L Fackler. Applied Computational
      Economics and Finance, MIT Press, 2002.
  
      """
      if n % 2 == 0:
          print("WARNING qnwsimp: n must be an odd integer. Increasing by 1")
          n += 1
  
      nodes = np.linspace(a, b, n)
      dx = nodes[1] - nodes[0]
      weights = np.kron(np.ones((n+1) // 2), np.array([2.0, 4.0]))
      weights = weights[:n]
      weights[0] = weights[-1] = 1
      weights = (dx / 3.0) * weights
  
      return nodes, weights