How to resolve zero division error for the following code which uses scipy?

Here's my python code:

import numpy as np
from scipy.integrate import quad

a = 0.15

def A(z):
    return -a * z**2

def integrand(x, z, B):
    return np.sqrt(-(2/x)*(3*x*A(x).derivative(n=2) - 3*x*(A(x).derivative(n=1))**2 + 6*A(x).derivative(n=1) + 2*B**4*x**3 + 2*B**2*x))

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

phi_values = [Phi(z, 0) for z in range(101)]

I am getting a "ZeroDivisionError: float division by zero"

How to resolve this issue?

I tried using sympy but to no avail.

Fixing your other breakage re. derivatives with some wild guessing, and adding an epsilon as your lower integration bound,

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]

But the saner thing to do is analytically divide out 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))

which yields the same results.