code: def hessian(self, x, y, kwargs, k=None):
          
          x = np.array(x, dtype=float)
          y = np.array(y, dtype=float)
          if self._foreground_shear:
              # needs to be computed numerically due to non-linear effects
              f_xx, f_xy, f_yx, f_yy = self.hessian_differential(x, y, kwargs, k=k)
          else:
              bool_list = self._bool_list(k)
              x_ = x
              y_ = y
              f_xx, f_yy, f_xy = np.zeros_like(x_), np.zeros_like(x_), np.zeros_like(x_)
              for i, func in enumerate(self.func_list):
                  if bool_list[i] is True:
                      f_xx_i, f_yy_i, f_xy_i = func.hessian(x_, y_, **kwargs[i])
                      f_xx += f_xx_i
                      f_yy += f_yy_i
                      f_xy += f_xy_i
          f_yx = f_xy
          return f_xx, f_xy, f_yx, f_yy
docs: hessian matrix
          :param x: x-position (preferentially arcsec)
          :type x: numpy array
          :param y: y-position (preferentially arcsec)
          :type y: numpy array
          :param kwargs: list of keyword arguments of lens model parameters matching the lens model classes
          :param k: only evaluate the k-th lens model
          :return: f_xx, f_xy, f_yy components