Btukfyl

Hi I'm trying to sove the Schrödinger equation for a infinite square well with a barrior using solve_bvp but I don't understand how to do two things with solve_bvp. one is setting the range, the x axis I am trying to give it goes from -b to b but it is solving for 0 to b, I had previously tried solving with solve_ivp but I was running into the same problem. the other problem I have is how to set the boundary condition, I need to take into account the first and second derivative, but I'm a bit lost because I cannot understand the documentation very well. I was following an example found here, but it wasn't giving me the correct results so I'm trying to rewrite it. I'm also discussing the problem on another site but that one is a bit more general than what I'm trying to solve here.

this is the code that I'm trying to get to work

from scipy.integrate import solve_bvp

from pylab import *

import numpy as np



## dimentions         _

## |       :       |   | Vmax

## |    ___:___    |   |        _

## |___|       |___|  _| Vmin   _|- Vb

##-B  -A   0   A   B*Alpha

## ______/|______/

##     Po       Pf



Vmax = 55

Vmin = 0

Vb = 50

A = 1

B = 4

alpha = 1

Po = -A-B

Pf = A+B*alpha

psi_b = array([0,1])   # Wave function boundry [first dir, secound dir]



N = 1000 #steps

psi = np.zeros([N,2])

x = linspace(Po, Pf, N)    # x-axis





def V(z):

    '''

    #Potential function in the finite square well.

    '''

    val=

    for i in z:

        if -A <=i <= A:

            val.append( Vb)

        elif i<=Po:

            val.append( Vmax)

        elif i>=Pf:

            val.append(Vmax)

        else:

            val.append(Vmin)

    return val



def boundary_conditions(psi_Po, psi_pf):

    return (psi_b)



def SE(z, p, E):

    state0 = p[1]

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

    return np.array([state0, state1])



def Wave_function(energy,psi):

    y_ = np.zeros((2, x.size))

    psi = solve_bvp(SE(x, psi, y_), boundary_conditions, x, y_)

    print(len(psi.y))

    return psi.y



def pltPotentail():

    plot(x,V(x))    



def main():

    en = linspace(0, Vb, 10)   # vector of energies where we look for the stable states

    figure(1)

    pltPotentail()

    Wave_function(0,psi)

    show()

if __name__ == "__main__":

    main()