code: def roots(self):
          
          if self.size() == 1:
              return np.array([])
  
          elif self.size() <= 100:
              ak = self.coefficients()
              v = np.zeros_like(ak[:-1])
              v[1] = 0.5
              C1 = linalg.toeplitz(v)
              C2 = np.zeros_like(C1)
              C1[0,1] = 1.
              C2[-1,:] = ak[:-1]
              C = C1 - .5/ak[-1] * C2
              eigenvalues = linalg.eigvals(C)
              roots = [eig.real for eig in eigenvalues
                      if np.allclose(eig.imag,0,atol=1e-10)
                          and np.abs(eig.real) <=1]
              scaled_roots = self._ui_to_ab(np.array(roots))
              return scaled_roots
          else:
              # divide at a close-to-zero split-point
              split_point = self._ui_to_ab(0.0123456789)
              return np.concatenate(
                  (self.restrict([self._domain[0],split_point]).roots(),
                   self.restrict([split_point,self._domain[1]]).roots())
              )
docs: Utilises Boyd's O(n^2) recursive subdivision algorithm. The chebfun
          is recursively subsampled until it is successfully represented to
          machine precision by a sequence of piecewise interpolants of degree
          100 or less. A colleague matrix eigenvalue solve is then applied to
          each of these pieces and the results are concatenated.
          See:
          J. P. Boyd, Computing zeros on a real interval through Chebyshev
          expansion and polynomial rootfinding, SIAM J. Numer. Anal., 40
          (2002), pp. 1666–1682.