docs: def interpolation_points(self, N):
          """
          N Chebyshev points in [-1, 1], boundaries included
          """
          if N == 1:
              return np.array([0.])
          return np.cos(np.arange(N)*np.pi/(N-1))