用一些疯狂的猜测来修正你的其他断裂re.导数，并添加一个ε作为你的积分下限，

import numpy as np
from scipy.integrate import quad

a = 0.15

# def A(z: float) -> float:
#     return -a * z**2
A = np.poly1d((-a, 0, 0))
Ad1 = A.deriv(m=1)
Ad2 = A.deriv(m=2)

def integrand(x: float, z: float, B: float) -> float:
    return np.sqrt(
        -(2/x)*(
            3*x*Ad2(x)
            - 3*x*Ad1(x)**2
            + 6*Ad1(x)
            + 2*B**4*x**3
            + 2*B**2*x
        )
    )

def Phi(z: float, B: float):
    epsilon = 1e-6
    result, _ = quad(func=integrand, a=epsilon, b=z, args=(z, B))
    return result

phi_values = [Phi(z, 0) for z in range(101)]
print(np.array(phi_values))

[-2.32379001e-06  2.36195636e+00  4.94106004e+00  7.90800675e+00
  1.13771987e+01  1.54207025e+01  2.00839857e+01  2.53965776e+01
...
  3.40309939e+03  3.47405018e+03  3.54573543e+03  3.61815514e+03
  3.69130933e+03]

但更明智的做法是分析性地划分x：

import numpy as np
from numpy.polynomial import Polynomial
from scipy.integrate import quad

a = 0.15

# def A(z: float) -> float:
#     return -a * z**2
A = Polynomial((0, 0, -a))
Ad1 = A.deriv(m=1)
Ad2 = A.deriv(m=2)

# Divide by x prior to integrand
Ad1ox = Ad1 // Polynomial((0, 1))

def integrand(x: float, z: float, B: float) -> float:
    return np.sqrt(
        -2*(
            3*Ad2(x)
            - 3*Ad1(x)**2
            + 6*Ad1ox(x)  # 6/x*Ad1(x)
            + 2*B**4*x**2
            + 2*B**2
        )
    )

def Phi(z: float, B: float):
    result, _ = quad(func=integrand, a=0, b=z, args=(z, B))
    return result

phi_values = [Phi(z, 0) for z in range(101)]
print(np.array(phi_values))

这产生相同的结果。