code: def jac(self, xy=None):
          
          if xy is not None:
              self.set_xy(xy)
  
          ux = numpy.dot(self.ax, self.xy_list[: len(self.ax)])
          vx = numpy.dot(self.bx, self.xy_list[: len(self.bx)])
          uy = numpy.dot(self.ay, self.xy_list[: len(self.ay)])
          vy = numpy.dot(self.by, self.xy_list[: len(self.by)])
  
          ux_dx = numpy.dot(self.ax, self.dx_list[: len(self.ax)])
          vx_dx = numpy.dot(self.bx, self.dx_list[: len(self.bx)])
          uy_dx = numpy.dot(self.ay, self.dx_list[: len(self.ay)])
          vy_dx = numpy.dot(self.by, self.dx_list[: len(self.by)])
  
          ux_dy = numpy.dot(self.ax, self.dy_list[: len(self.ax)])
          vx_dy = numpy.dot(self.bx, self.dy_list[: len(self.bx)])
          uy_dy = numpy.dot(self.ay, self.dy_list[: len(self.ay)])
          vy_dy = numpy.dot(self.by, self.dy_list[: len(self.by)])
  
          jac = numpy.array(
              [
                  [
                      (ux_dx * vx - vx_dx * ux) / vx ** 2,
                      (ux_dy * vx - vx_dy * ux) / vx ** 2,
                  ],
                  [
                      (uy_dx * vy - vy_dx * uy) / vy ** 2,
                      (uy_dy * vy - vy_dy * uy) / vy ** 2,
                  ],
              ]
          )
          return jac
docs: Get the Jacobian at (x, y).