code: def hessian(self, x, y, theta_E, center_x=0, center_y=0):
          
          x_shift = x - center_x
          y_shift = y - center_y
          R = (x_shift*x_shift + y_shift*y_shift)**(3./2)
          if isinstance(R, int) or isinstance(R, float):
              prefac = theta_E / max(0.000001, R)
          else:
              prefac = np.empty_like(R)
              r = R[R>0]  #in the SIS regime
              prefac[R==0] = 0.
              prefac[R>0] = theta_E / r
  
          f_xx = y_shift*y_shift * prefac
          f_yy = x_shift*x_shift * prefac
          f_xy = -x_shift*y_shift * prefac
          return f_xx, f_yy, f_xy
docs: returns Hessian matrix of function d^2f/dx^2, d^f/dy^2, d^2/dxdy