docs: def _qnwbeta1(n, a=1.0, b=1.0):
      """
      Computes nodes and weights for quadrature on the beta distribution.
      Default is a=b=1 which is just a uniform distribution
  
      NOTE: For now I am just following compecon; would be much better to
      find a different way since I don't know what they are doing.
  
      Parameters
      ----------
      n : scalar : int
          The number of quadrature points
  
      a : scalar : float, optional(default=1)
          First Beta distribution parameter
  
      b : scalar : float, optional(default=1)
          Second Beta distribution parameter
  
      Returns
      -------
      nodes : np.ndarray(dtype=float, ndim=1)
          The quadrature points
  
      weights : np.ndarray(dtype=float, ndim=1)
          The quadrature weights that correspond to nodes
  
      Notes
      -----
      Based of original function ``_qnwbeta1`` in CompEcon toolbox by
      Miranda and Fackler
  
      References
      ----------
      Miranda, Mario J, and Paul L Fackler. Applied Computational
      Economics and Finance, MIT Press, 2002.
  
      """
      # We subtract one and write a + 1 where we actually want a, and a
      # where we want a - 1
      a = a - 1
      b = b - 1
  
      maxiter = 25
  
      # Allocate empty space
      nodes = np.zeros(n)
      weights = np.zeros(n)
  
      # Find "reasonable" starting values.  Why these numbers?
      for i in range(n):
          if i == 0:
              an = a/n
              bn = b/n
              r1 = (1+a) * (2.78/(4+n*n) + .768*an/n)
              r2 = 1 + 1.48*an + .96*bn + .452*an*an + .83*an*bn
              z = 1 - r1/r2
          elif i == 1:
              r1 = (4.1+a) / ((1+a)*(1+0.156*a))
              r2 = 1 + 0.06 * (n-8) * (1+0.12*a)/n
              r3 = 1 + 0.012*b * (1+0.25*abs(a))/n
              z = z - (1-z) * r1 * r2 * r3
          elif i == 2:
              r1 = (1.67+0.28*a)/(1+0.37*a)
              r2 = 1+0.22*(n-8)/n
              r3 = 1+8*b/((6.28+b)*n*n)
              z = z-(nodes[0]-z)*r1*r2*r3
          elif i == n - 2:
              r1 = (1+0.235*b)/(0.766+0.119*b)
              r2 = 1/(1+0.639*(n-4)/(1+0.71*(n-4)))
              r3 = 1/(1+20*a/((7.5+a)*n*n))
              z = z+(z-nodes[-4])*r1*r2*r3
          elif i == n - 1:
              r1 = (1+0.37*b) / (1.67+0.28*b)
              r2 = 1 / (1+0.22*(n-8)/n)
              r3 = 1 / (1+8*a/((6.28+a)*n*n))
              z = z+(z-nodes[-3])*r1*r2*r3
          else:
              z = 3*nodes[i-1] - 3*nodes[i-2] + nodes[i-3]
  
          ab = a+b
  
          # Root finding
          its = 0
          z1 = -100
          while abs(z - z1) > 1e-10 and its < maxiter:
              temp = 2 + ab
              p1 = (a-b + temp*z)/2
              p2 = 1
  
              for j in range(2, n+1):
                  p3 = p2
                  p2 = p1
                  temp = 2*j + ab
                  aa = 2*j * (j+ab)*(temp-2)
                  bb = (temp-1) * (a*a - b*b + temp*(temp-2) * z)
                  c = 2 * (j - 1 + a) * (j - 1 + b) * temp
                  p1 = (bb*p2 - c*p3)/aa
  
              pp = (n*(a-b-temp*z) * p1 + 2*(n+a)*(n+b)*p2)/(temp*(1 - z*z))
              z1 = z
              z = z1 - p1/pp
  
              if abs(z - z1) < 1e-12:
                  break
  
              its += 1
  
          if its == maxiter:
              raise ValueError("Max Iteration reached.  Failed to converge")
  
          nodes[i] = z
          weights[i] = temp/(pp*p2)
  
      nodes = (1-nodes)/2
      weights = weights * math.exp(gammaln(a+n) + gammaln(b+n)
                                   - gammaln(n+1) - gammaln(n+ab+1))
      weights = weights / (2*math.exp(gammaln(a+1) + gammaln(b+1)
                           - gammaln(ab+2)))
  
      return nodes, weights