x_values = arange(0.01, 0.1, 0.0001)
delta = []
for x in x_values:

    def fun(E):
        distance = np.sqrt(E * E + x * x)
return np.tanh(1477.92 * distance) / distance
integral = integrate.quad(fun, 0, 1)
delta_val = fsolve(integral, 1e-23) - a
delta.append(delta_val)
pl.plot(delta, x_values)
Find root of numerical integration