import numpy as np

import matplotlib.pyplot as plt

import scipy.optimize as opt

import scipy.integrate as inte

# Variables.

boltzmann_const = 1.38e-23

planck_const = 6.62e-34

hbar = planck_const / ( 2 * np.pi )

transition_temp = 9.2

gap_energy_at_zero_kelvin = 3.528 / ( 2 * transition_temp * boltzmann_const )

debye_freq = ( 296 * boltzmann_const ) / hbar


# For subtracting from root_of_integral

a_const = np.log( ( 1.13 * hbar * debye_freq ) / ( boltzmann_const * transition_temp) )

# For simplifying function f.

b_const = ( hbar * debye_freq ) / ( 2 * boltzmann_const)


x_values = np.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 = inte.quad( fun, 0, 1 )

    delta_val = opt.fsolve(lambda x:integral ,1e-23 ) - a_const

    delta.append( delta_val )

    plt.plot( delta,x_values )




ValueError: x and y must have same first dimension, but have shapes (1, 1) and (901,)