title: Find root of numerical integration
body: <p>I am trying to reproduce this Mathematica program in Python:</p>
  
  <p><img src="https://i.stack.imgur.com/0Km6p.png" alt="Mathematica program"></p>
  
  <p>It finds the roots of an numerical integration, and forms a plot of these values. However, I cannot get my attempt to run.</p>
  
  <p>Current attempt:</p>
  
  <p>from scipy.integrate import quad
  from scipy import integrate
  from scipy.optimize import fsolve
  import pylab as pl
  import numpy as np</p>
  
  <pre><code># 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)
  
  
  def f( coherence_length, temp ):
      # Defines the equation whose integral will have its roots found. Epsilon = coherence length. Delta = Gap energy.
  
      squareRoot = np.sqrt( coherence_length*coherence_length + gap_energy*gap_energy )
      return np.tanh( ( ( b_const / temp ) * squareRoot ) / squareRoot )
  
  
  def integrate( coherence_length, temp ):
      # Integrates equation f with respect to E, between 0 and 1. 
  
      return integrate.quad( f, 0, 1, args = ( temp, ) )[0]
  
  
  def root_of_integral( temp ):
      # Finds the roots of the integral with a guess of 0.01.
  
     return fsolve( integrate, 0.01, args = ( temp, ) )
  
  
  def gap_energy_values( temp ):
      # Subtracts a_const from each root found, to obtain the gap_energy_values.
  
      return root_of_integral( temp ) - a_const
  </code></pre>