How To Solve A Simple Boundary Value Problem For Tise On Python

I am trying to solve the TISE for an infinite potential well V=0 on the interval [0,L]. The exercise gives us that the value of the wavefunction and its derivative at 0 is 0,1 resp

Solution 1:

For all the other solvers in scipy the argument order x,y, and even in odeint one can use this order by giving the option tfirst=True. Thus change to

defeqn(x, y, energy):              #array of first order ODE's     
    y0, y1 = y
    y2 = -2*m_el*energy*y0/hbar**2return [y1,y2]

For the BVP solver you have to think of the energy parameter as an extra state component with zero derivative, thus adding a third slot in the boundary conditions. Scipy's solve_bvp allows to keep it as parameter, so that you get 3 slots in the boundary conditions, allowing to fix the first derivative at x=0 to select one non-trivial solution from the eigenspace.

defbc(y0, yL, E):
    return[ y0[0], y0[1]-1, yL[0] ]

Next construct an initial state that is close to the suspected ground state and call the solver

x0 = np.linspace(0,L_bohr,6);
y0 = [ x0*(1-x0/L_bohr), 1-2*x0/L_bohr ]
E0 = 134*e_el

sol = solve_bvp(eqn, bc, x0, y0, p=[E0])
print(sol.message, "  E=", sol.p[0]/e_el," eV")

and then produce the plot

x = np.linspace(0,L_bohr,1000)
plt.plot(x/L_bohr, sol.sol(x)[0]/L_bohr,'-+', ms=1)
plt.grid()

The algorithm converged to the desired accuracy. E= 134.29310361903723 eV