code: def _2ndDerivInt(x,y,z,a2,b2,c2,n,i,j):
      
      def integrand(tau):
          if i!=j:
              return _FracInt(x,y,z,a2,b2,c2,tau,n-1)*n*\
                  (1.+(-1.-2.*x/(tau+a2))*(i==0 or j==0))*\
                  (1.+(-1.-2.*y/(tau+b2))*(i==1 or j==1))*\
                  (1.+(-1.-2.*z/(tau+c2))*(i==2 or j==2))
          else:
              var2 = x**2*(i==0) + y**2*(i==1) + z**2*(i==2)
              coef2 = a2*(i==0) + b2*(i==1) + c2*(i==2)
              return _FracInt(x,y,z,a2,b2,c2,tau,n-1)*n*(4.*var2)/(tau+coef2)**2 +\
                  _FracInt(x,y,z,a2,b2,c2,tau,n)*(-2./(tau+coef2))
      return integrate.quad(integrand, lowerlim(x**2,y**2,z**2,a2,b2,c2), np.inf)[0]
docs: Integral involved in second derivatives d^\Phi/(dx_i dx_j)
      integrate
          1/A B^(n-1) (-2 x_i/(tau+a_i)) (-2 x_j/(tau+a_j))
      when i /= j or
          1/A [ B^(n-1) 4n x_i^2 / (a_i+t)^2 + B^n -(-2/(a_i+t)) ]
      when i == j where
      A = sqrt((tau+a)(tau+b)(tau+c)) and B = (1-x^2/(tau+a)-y^2/(tau+b)-z^2/(tau+c))
      from lambda to infty with respect to tau
      The lower limit lambda is given by lowerlim function.
      This is a second derivative of _potInt.