Agfdhyk

I've been trying to use either solve_ivp or solve_bvp to solve a problem I've been having but I'm not making any progress. I think the code I have here will work but I cannot get the range to be correct. for some reason I cannot understand the range is always going from 0 to x and not from -x to x can someone help me fix this part for solve_ivp?

here is the code reduced to a min

from pylab import *

from scipy.integrate import solve_ivp

from scipy.optimize import brentq

import numpy as np

import itertools



a=1

B=4

L= B+a

Vmax= 50

Vpot = False



N = 1000                  # number of points to take

psi = np.zeros([N,2])     # Wave function values and its derivative (psi and psi')

psi0 = array([0,1])   # Wave function initial states

Vo = 50

E = 0.0                   # global variable Energy  needed for Sch.Eq, changed in function "Wave function"

b = L                     # point outside of well where we need to check if the function diverges

x = linspace(-B-a, L, N)    # x-axis

def V(x):

    '''

    #Potential function in the finite square well.

    '''

    if -a <=x <=a:

        val = Vo

    elif x<=-a-B:

        val = Vmax

    elif x>=L:

        val = Vmax

    else:

        val = 0

    return val



def SE(z, p):

    state0 = p[1]

    state1 = 1.0*(V(z) - E)*p[0]

    return array([state0, state1])



def Wave_function(energy):

    global E

    E = energy

    #        odeint(func, y0, t)

    #     solve_ivp(fun, t_span, y0)

    psi = solve_ivp(SE, [-B-a, L], psi0, max_step = ((B+a+L)/(N)))

    return psi.y



def main():

    # main program        

    f2 = figure(2)

    plot(x, Wave_function(9.8)[0][:1000])

    grid()

    f2.show()



if __name__ == "__main__":

    main()

this is what the code gives me in the end. the right side is okay but the left side is wrong. I depend on both sides working, not for visuals.

edit:
for charity this is what the potential function should look like:

the final graph should look similar to this:

Not sure if it helps but it may gives you a hint.
It is not that solve_ivp doesn't work for -x to 0, but your function V may be wrong. I noticed that the wave begins to appear after that V decreases from Vmax to 0.

This code :

%matplotlib inline

from pylab import *

from scipy.integrate import solve_ivp

from scipy.optimize import brentq

import numpy as np

import itertools



a=1.

B=4.

L= B+a

Vmax= 50.



N = 10000

E = 9.8



def V(x):

    if -L <= x <= -B:

        return Vmax

    else:

        return 0



def SE(z, p):

    state0 = p[1]

    state1 = (V(z) - E)*p[0]

    return array([state0, state1])



def Wave_function():

    return solve_ivp(SE, [-L, L], [0., 1.], max_step=2*L/N)



result = Wave_function()

plot(result.t, result.y[0], color='tab:blue')

gives you the "expected" output :

Your code seems Okay in general. However, given your figure for the potential energy, the value for Vo should be *Vo = 10. In addition, in your main function your are only plotting the wave function as the solution of the Schrodinger Equation. Bellow, is what I am proposing you as a possible solution to your problem assuming that I properly understood your concern:

import matplotlib.pyplot as plt

from scipy.integrate import solve_ivp

import numpy as np



a=1

B=4

L= B+a

Vmax= 50





N = 1000                  # number of points to take

psi = np.zeros([N,2])     # Wave function values and its derivative (psi and psi')

psi0 = np.array([0,1])   # Wave function initial states

Vo = 10     # Not 50, in order to conform your figure of the potential energy

E = 0.0         # global variable Energy  needed for Sch.Eq, changed in    

                # function "Wave function"

b = L           # point outside of well where we need to check if the function diverges

x = np.linspace(-L, L, N) # linspace(-B-a, L, N)    # x-axis





def V(x):

    '''

    Potential function in the finite square well.

    '''

    if -a <=x <=a:

        val = Vo

    elif x<= -L:  # -a-B:

        val = Vmax

    elif x>=L:

        val = Vmax

    else:

        val = 0

    return val



def SE(z, p):

    state0 = p[1]

    state1 = 1.0*(V(z) - E)*p[0]

    return array([state0, state1])



def Wave_function(energy):

    global E

    E = energy

    psi = solve_ivp(SE, [-B-a, L], psi0, max_step = ((B+a+L)/(N)))

    return psi.y



def main():

    # main program       

    plt.figure()

    plt.subplot(121)

    plt.plot(x, Wave_function(9.8)[0][:1000])

    plt.grid()

    plt.title("Wave function")

    plt.xlabel(r"$ x $")

    plt.ylabel(r"$psi(x)$")

    plt.subplot(122)



    potential = np.vectorize(V) # Make the function 'V(x)' to also work on array



    pot = potential(x) # Potential ernergy in the defined domain of 'x'

    t = [-L, -a, a, L] # the singular value of x

    y = potential(t)   # the potential energy at thos singular value of 'x'

    # But to conform your figure I'll just do y = 0 * y

    plt.plot(x, pot, t, 0*y, 'ko')

    plt.title("Potential Energy")

    plt.xlabel(r"$ x $")

    plt.ylabel(r"$V(x)$")

    plt.show()



if __name__ == "__main__":

    main()

The output figure is the following: