code: def plot(self, xmin=-1, xmax=1,
               resolution_outside=20, resolution_inside=200):
          
          if xmin > self.L or xmax < self.R:
              raise ValueError('xmin=%g > L=%g or xmax=%g < R=%g is meaningless for plot' % (xmin, self.L, xmax, self.R))
  
          if self.eps_L == 0 and self.eps_R == 0:
              return ([xmin, self.L, self.L, self.R, self.R, xmax],
                      [0, 0, 1, 1, 0, 0])
          else:
              n = float(resolution_inside)/(0.5*(self.eps_L + self.eps_R))
              x = np.concatenate((
                  np.linspace(xmin, self.L-self.eps_L, resolution_outside+1),
                  np.linspace(self.L-self.eps_L, self.R+self.eps_R,
                              resolution_inside+1),
                  np.linspace(self.R+self.eps_R, xmax, resolution_outside+1)))
              y = self(x)
              return x, y
docs: Return arrays x, y for plotting IndicatorFunction
          on [`xmin`, `xmax`]. For the exact discontinuous
          indicator function, we typically have
          ``x = [xmin, L, L, R, R, xmax]; y = [0, 0, 1, 1, 0, 0]``,
          while for the smoothed version, the densities of
          coordinates in the ``x`` array is computed on basis of the
          `eps` parameter with `resolution_outside` plotting intervals
          outside the smoothed regions and  `resolution_inside` intervals
          inside the smoothed regions.