code: def hessian(self, theta_x, theta_y, kwargs_lens, k=None, diff=0.00000001):
          
  
          alpha_ra, alpha_dec = self.alpha(theta_x, theta_y, kwargs_lens)
  
          alpha_ra_dx, alpha_dec_dx = self.alpha(theta_x + diff, theta_y, kwargs_lens)
          alpha_ra_dy, alpha_dec_dy = self.alpha(theta_x, theta_y + diff, kwargs_lens)
  
          dalpha_rara = (alpha_ra_dx - alpha_ra)/diff
          dalpha_radec = (alpha_ra_dy - alpha_ra)/diff
          dalpha_decra = (alpha_dec_dx - alpha_dec)/diff
          dalpha_decdec = (alpha_dec_dy - alpha_dec)/diff
  
          f_xx = dalpha_rara
          f_yy = dalpha_decdec
          f_xy = dalpha_radec
          f_yx = dalpha_decra
          return f_xx, f_xy, f_yx, f_yy
docs: computes the hessian components f_xx, f_yy, f_xy from f_x and f_y with numerical differentiation
  
          :param theta_x: x-position (preferentially arcsec)
          :type theta_x: numpy array
          :param theta_y: y-position (preferentially arcsec)
          :type theta_y: numpy array
          :param kwargs_lens: list of keyword arguments of lens model parameters matching the lens model classes
          :param diff: numerical differential step (float)
          :return: f_xx, f_xy, f_yx, f_yy