[Python] Solving TISE using solve_ivp and then a root finder

I am solving the TISE for x between 0 and L. I am given the boundary conditions that at x=0, phi(0)=0 and dphi(0)/dx=1. I also know that at x=L, phi(L)=0. I have used solve_ivp and a suggested initial energy of 134.0 which creates a correct graph.

from scipy.constants import c, m_e, e
me_ev = m_e*c**2/e # mass m_e in [eV]

l_bohr = 5.2918e-11 # Bohr radius [m]
convert_length = 5067730.179 # convert [m] to [1/eV] in natural units
l_ev = l_bohr * convert_length

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from scipy.optimize import root_scalar

# TISE
def eqn (x, y, energy):
me_ev = m_e*c**2/e
f_0 = y[1] # dphi/dx
f_1 = -2 * me_ev * energy * y[0] # d2phi/dx2, y[0] = phi
return np.array([f_0, f_1])

# using solve_ivp and the initial coonditions phi(0)=0 and dphi(0)/dx=1 to find phi(L)
def solve(energy, func):
#constants
l_bohr = 5.2918e-11 # Bohr radius [m]
convert_length = 5067730.179 # convert [m] to [1/eV] in natural units
l_ev = l_bohr * convert_length

# x goes form [0,L]
x = np.linspace(0., l_ev, 1000)

# initial conditions
phi0 = 0.0
dphi0dx = 1.0
init = np.array([phi0, dphi0dx])

answer = solve_ivp(func, (0., l_ev), init, t_eval=x, args=(energy,))
phiL = answer.y[0][-1]
return phiL


This is where the code starts to have problems

# calculating ground state energy:
# use fsolve to put different energies into solve(energy, func) such that phiL is minimized:
ground_energy = fsolve(solve, 134.0, args=(eqn,))
print(ground_energy)


When I do this I get the error: setting an array element with a sequence. The requested array has an inhomogeneous shape after 1 dimensions. The detected shape was (2,) + inhomogeneous part.

Does anybody know how to fix this?

The function solve returns phi(L) so that I can optimize the energy such that phi(L)=0. This is where I run into problems. I try use fsolve.

Continue reading...