# Numerically solving the Schrodinger Equation

This python notebook will show a simple method to find solutions to the Schrodinger equation using what is called the shooting method. Recall that the Schrodinger equation can be written as:

$$ -\frac{\hbar^2}{2 m} \frac{\partial^2 \psi}{\partial x^2} +(V(x)- E) \psi =0 $$

where $V(x)$ is the potential and both the energy $E$ and the wavefunction $\psi$ are unknown.
The fact that we do not know the eigenvalue $E$ means that we have to use trial and error to guess a value for $E$ and the numerically solve Schrodinger's equation to find $\psi(x)$. We will know if our solution is correct if the wavefunction approaches 0 for large $x$.

To proceed we first change our units so that $\hbar^2/2 m=1$ and then add in a new variable $\phi$ which represents the derivative of $\psi$ to give us the set of first order differential equations:

$$ \begin{align}
    \frac{d \psi}{d x} &= \phi(x) \\
    \frac{d \phi}{d x} &= (V(x)-E) \psi(x)
    \end{align} $$

Recall that we did the same trick in the Cyclotron Motion notebook to convert the 2nd order differential equations into a set of first order differential equations. This means that we can use the inbuilt functions such as 'odeint' to automatically solve the equations for us.

So first as usual we import the required libraries


```python
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt

from scipy.integrate import odeint
```

Next we define the range over which we want to find the solutions. For the example below this is taken to be between $-20< x <20$ and we split it up into 2000 points.


```python
start_position=-20.0
end_position=20   #set this to be the correct end point.
space_steps=2000
step=(end_position-start_position)/(space_steps-1)
space=np.linspace(start_position,end_position,space_steps)
```

Next we define our potential function. This can be an arbitrary function but here I have taken it to be a square well with a negative depth and lying symmetrically about $x=0$.


```python
def potential(x,depth,width):  #square well potential
    pot=0.0
    if abs(x)<=width/2:
        pot=depth
    else:
        pot=0.0
    return pot
```

Now we define the derivatives functions for the equations above. These are the equations that we want to solve.


```python
def shoot(yin,position,depth,width,energy):
    return np.array([yin[1],(potential(position,depth,width)-energy)*yin[0]])
```

Now we still don't know the initial values for $\psi(-20)$ and $\phi(-20)$ since we don't know the energy eigenvalue. However the Schrodinger equation is linear which means that we are free to set $\psi(-20)=1$ since if $\psi(x)$ is a solution so is $\alpha \psi(x)$ for any complex number $\alpha$. We also know that $\psi(x)$ decays exponentially if the energy is less than the potential and that the rate of decay is given by $V(x)-E$. Thus we can approximate the derivative of $psi(x)$ (i.e. $\phi(x)$) as $\phi(x) \approx (V(x)-E) \psi(x)$.

We then guess an energy level $E_g$ which means that the initially
$\phi(-20)=(V(-20)-E_g)$. Since for the square well potential $V(-20)=0$ this just becomes $\phi(-20)=-E_g$.

The Shooting method thus involves guessing an initial value of the energy and integrating the equations to see if it that value is correct. To see this in action let us use a potential well with a depth of -5 and a width of 30. This ensures that there are plenty of solutions to the Schrodinger equation. Our first guess will be that E=-4.9896. 


```python
energy_guess=-4.9896  #this is our initial guess.
start=np.array([1.0,-energy_guess])

xout=odeint(shoot,start,space,args=(-5.0,30,energy_guess))
xt=np.transpose(xout)
psi=xt[0]
psiprime=xt[1]

plt.figure()
plt.plot(space,psi,linewidth=3,label='wavefunction')
plt.plot(space,psiprime,linewidth=3,label='derivative')
plt.legend()
plt.show()
```

Notice that the solution is clearly wrong since it heads towards $-\infty$ for large $x$. So
we can amend our guess of the energy for something larger and try again.


```python
energy_guess=-4.9897  #this is our second guess.
start=np.array([1.0,-energy_guess])

xout=odeint(shoot,start,space,args=(-5.0,30,energy_guess))
xt=np.transpose(xout)
psi=xt[0]
psiprime=xt[1]

plt.figure()
plt.plot(space,psi,linewidth=3,label='wavefunction')
plt.plot(space,psiprime,linewidth=3,label='derivative')
plt.legend()
plt.show()

```

Now the solution is still wrong but it heads towards $+\infty$ for large $x$. We can thus say that the correct value for the energy must lie between our two solutions and all we have to do is find it. Trial and error leads to the solution below.


```python
energy_guess=-4.989659756171312  #this is ground state
start=np.array([1.0,-energy_guess])

xout=odeint(shoot,start,space,args=(-5.0,30,energy_guess))
xt=np.transpose(xout)
psi=xt[0]
psiprime=xt[1]

plt.figure()
plt.plot(space,psi,linewidth=3,label='wavefunction')
plt.plot(space,psiprime,linewidth=3,label='derivative')
plt.legend()
plt.show()
```

Note that as expected the ground state wavefunction has a single maximum and goes to zero as the distance from the well increases. If we wanted to we could normalise the wavefunction so that the area under the curve was one but for now we will look at other solutions. A bit more trial and error leads to the next guess for the energy of the first higher order state.


```python
energy_guess=-4.95864151237514  #this is first higher order state
start=np.array([1.0,-energy_guess])

xout=odeint(shoot,start,space,args=(-5.0,30,energy_guess))
xt=np.transpose(xout)
psi=xt[0]
psiprime=xt[1]
plt.figure()
plt.plot(space,psi,linewidth=3,label='wavefunction')
plt.plot(space,psiprime,linewidth=3,label='derivative')
plt.legend()
plt.show()
```

It very quickly gets complicated to continually bisect the interval to find new
solutions. However it is fairly easy to write a program that does that for us. There are lots of different ways of doing so since it is the same as finding the zero of a function but this is perhaps the simplest. We start with two values $E_1$ and $E_2$ which have the property that
the solution to the equation goes to $\pm \infty$ for the two values and then we look at the
sign of the solution at the midpoint.

If the sign of the solution at the midpoint is the same value as $E_1$ we replace $E_1$ by
$(E_1+E_2)/2$ and start again. Otherwise we replace $E_2$ by the midpoint and start again. 
Each time we do this the size of the interval decreases by a factor of 2 so we can continue doing it for as long as we like. In practice we stop when we run out of numerical precision.



```python
def find(start,stop): #make sure that the solution goes in different directions for the two values of start and stop
    out1=odeint(shoot,np.array([1.0,start]),space,args=(-5.0,30,start))[-1][0]
    #out1 is the final value of our trial solution for the guess Energy=start
    # the element [-1][0] picks the last element in the array.
    out2=odeint(shoot,np.array([1.0,stop]),space,args=(-5.0,30,stop))[-1][0]
    #out2 is the final value of our trial solution for the guess Energy=stop.
    guess=(stop+start)/2.0
    #guess is the midpoint of the interval.
    out3=odeint(shoot,np.array([1.0,guess]),space,args=(-5.0,30,guess))[-1][0]
    #out3 is the final value of the trial solution for the guess Energy=guess.
    print(out1,out2,out3) #prints out the three values. Check that two of them have different signs.
    count=1
    while count<50:
        if np.sign(out1)==np.sign(out3):
            start=guess
        else:
            stop=guess
        out1=odeint(shoot,np.array([1.0,start]),space,args=(-5.0,30,start))[-1][0]
        out2=odeint(shoot,np.array([1.0,stop]),space,args=(-5.0,30,stop))[-1][0]
        guess=(stop+start)/2.0
        out3=odeint(shoot,np.array([1.0,guess]),space,args=(-5.0,30,guess))[-1][0]
        print(out1,out2,out3)
        count=count+1
    return guess
```

Let us see it in action


```python
find(-4.90,-4.91)
```

    3017157085.367703 -1488721157.6751328 907699261.0260235
    907699261.0260235 -1488721157.6751328 -261401555.8847637
    907699261.0260235 -261401555.8847637 331371480.6051831
    331371480.6051831 -261401555.8847637 36928602.6327268
    36928602.6327268 -261401555.8847637 -111764989.27813512
    36928602.6327268 -111764989.27813512 -37298344.358756356
    36928602.6327268 -37298344.358756356 -154677.68537378384
    36928602.6327268 -154677.68537378384 18394540.219054118
    18394540.219054118 -154677.68537378384 9121623.380339025
    9121623.380339025 -154677.68537378384 4484022.385920381
    4484022.385920381 -154677.68537378384 2164787.3561655935
    2164787.3561655935 -154677.68537378384 1004925.6042926874
    1004925.6042926874 -154677.68537378384 424972.10375108715
    424972.10375108715 -154677.68537378384 135111.13238121656
    135111.13238121656 -154677.68537378384 -9770.428800546533
    135111.13238121656 -9770.428800546533 62821.2645162628
    62821.2645162628 -9770.428800546533 26474.387481207268
    26474.387481207268 -9770.428800546533 8541.710969385915
    8541.710969385915 -9770.428800546533 -672.443897177723
    8541.710969385915 -672.443897177723 3904.8691690621663
    3904.8691690621663 -672.443897177723 1737.0947151509438
    1737.0947151509438 -672.443897177723 457.96044886326
    457.96044886326 -672.443897177723 -191.7039070536786
    457.96044886326 -191.7039070536786 -1.0829751041276638
    457.96044886326 -1.0829751041276638 408.9111314512371
    408.9111314512371 -1.0829751041276638 338.76562397636997
    338.76562397636997 -1.0829751041276638 370.21726566528537
    370.21726566528537 -1.0829751041276638 194.1052605534261
    194.1052605534261 -1.0829751041276638 102.73788347335487
    102.73788347335487 -1.0829751041276638 145.53026142253935
    145.53026142253935 -1.0829751041276638 126.99669441866197
    126.99669441866197 -1.0829751041276638 194.37546577396165
    194.37546577396165 -1.0829751041276638 61.788098655291165
    61.788098655291165 -1.0829751041276638 126.7005616086945
    126.7005616086945 -1.0829751041276638 170.39041778694593
    170.39041778694593 -1.0829751041276638 102.72423431724332
    102.72423431724332 -1.0829751041276638 149.3392670058012
    149.3392670058012 -1.0829751041276638 321.5967050125969
    321.5967050125969 -1.0829751041276638 168.24966070965803
    168.24966070965803 -1.0829751041276638 93.45030042691742
    93.45030042691742 -1.0829751041276638 278.85840376263917
    278.85840376263917 -1.0829751041276638 303.66334917688056
    303.66334917688056 -1.0829751041276638 104.22018498582312
    104.22018498582312 -1.0829751041276638 -1.0829751041276638
    104.22018498582312 -1.0829751041276638 -1.0829751041276638
    104.22018498582312 -1.0829751041276638 -1.0829751041276638
    104.22018498582312 -1.0829751041276638 -1.0829751041276638
    104.22018498582312 -1.0829751041276638 -1.0829751041276638
    104.22018498582312 -1.0829751041276638 -1.0829751041276638
    104.22018498582312 -1.0829751041276638 -1.0829751041276638





    -4.906952798962594




```python
energy_guess=-4.906952798962594  #this is 2nd higher order state
start=np.array([1.0,-energy_guess])

xout=odeint(shoot,start,space,args=(-5.0,30,energy_guess))
xt=np.transpose(xout)
psi=xt[0]
psiprime=xt[1]
plt.figure()
plt.plot(space,psi,linewidth=3,label='wavefunction')
plt.plot(space,psiprime,linewidth=3,label='derivative')
plt.legend()
plt.show()
```


```python
find(-4.80,-4.90)
```

    -4722982155.803635 3017157085.367703 3644517724.8830185
    -4722982155.803635 3644517724.8830185 -1986229942.4400694
    -1986229942.4400694 3644517724.8830185 664398818.885519
    -1986229942.4400694 664398818.885519 -737706816.3640819
    -737706816.3640819 664398818.885519 -52098561.277820334
    -52098561.277820334 664398818.885519 302820375.6614286
    -52098561.277820334 302820375.6614286 124459119.31465241
    -52098561.277820334 124459119.31465241 35946462.74768223
    -52098561.277820334 35946462.74768223 -8135558.151578826
    -8135558.151578826 35946462.74768223 13890764.072653009
    -8135558.151578826 13890764.072653009 2873808.451028525
    -8135558.151578826 2873808.451028525 -2631632.931605917
    -2631632.931605917 2873808.451028525 120813.50348482178
    -2631632.931605917 120813.50348482178 -1255499.5788772139
    -1255499.5788772139 120813.50348482178 -567478.2433713758
    -567478.2433713758 120813.50348482178 -223271.08395805047
    -223271.08395805047 120813.50348482178 -51322.461502839724
    -51322.461502839724 120813.50348482178 34940.31432688594
    -51322.461502839724 34940.31432688594 -8289.316584468526
    -8289.316584468526 34940.31432688594 13181.111317740812
    -8289.316584468526 13181.111317740812 2436.478239705172
    -8289.316584468526 2436.478239705172 -2842.9567906404345
    -2842.9567906404345 2436.478239705172 -195.00066379703512
    -195.00066379703512 2436.478239705172 1137.0215657811093
    -195.00066379703512 1137.0215657811093 417.8728060291689
    -195.00066379703512 417.8728060291689 79.18486656656538
    -195.00066379703512 79.18486656656538 -136.9772077683433
    -136.9772077683433 79.18486656656538 94.52989337020101
    -136.9772077683433 94.52989337020101 -52.6522230368047
    -52.6522230368047 94.52989337020101 226.26336513107776
    -52.6522230368047 226.26336513107776 15.282280136682667
    -52.6522230368047 15.282280136682667 189.10340411243791
    -52.6522230368047 189.10340411243791 -39.14502692750083
    -39.14502692750083 189.10340411243791 -35.35989478974898
    -35.35989478974898 189.10340411243791 -29.503347325512806
    -29.503347325512806 189.10340411243791 -3.6567437890452865
    -3.6567437890452865 189.10340411243791 -51.45818994412625
    -51.45818994412625 189.10340411243791 11.943022793374826
    -51.45818994412625 11.943022793374826 16.325568984438576
    -51.45818994412625 16.325568984438576 -58.55372925590135
    -58.55372925590135 16.325568984438576 -80.92186632553381
    -80.92186632553381 16.325568984438576 -33.33514981250881
    -33.33514981250881 16.325568984438576 -57.61579674889085
    -57.61579674889085 16.325568984438576 14.02191338249445
    -57.61579674889085 14.02191338249445 7.272817529658203
    -57.61579674889085 7.272817529658203 -49.624550862754404
    -49.624550862754404 7.272817529658203 16.67004679117914
    -49.624550862754404 16.67004679117914 -49.624550862754404
    -49.624550862754404 16.67004679117914 -49.624550862754404
    -49.624550862754404 16.67004679117914 -49.624550862754404





    -4.834606398105997




```python
energy_guess=-4.834606398105997  #which state is this?
start=np.array([1.0,-energy_guess])

xout=odeint(shoot,start,space,args=(-5.0,30,energy_guess))
xt=np.transpose(xout)
psi=xt[0]
psiprime=xt[1]
plt.figure()
plt.plot(space,psi,linewidth=3,label='wavefunction')
plt.plot(space,psiprime,linewidth=3,label='derivative')
plt.legend()
plt.show()
```

## Simple Harmonic Potential

In class we have also looked at a harmonic potential well defined by

$$ V(x) =\frac{x^2}{w^2} $$

where $w$ defines the steepness of the well. We can easily adapt our code to find the wavefunctions for this case as well.  


```python
def potential(x,width):  #harmonic potential well
    return x*x/(width*width)

def shoot(yin,position,width,energy):
    return np.array([yin[1],(potential(position,width)-energy)*yin[0]])
```


```python
def find_harmonic(start,stop,width): #make sure that the solution goes in different directions for the two values of start and stop
    out1=odeint(shoot,np.array([1.0,potential(-20,width)-start]),space,args=(width,start))[-1][0]
    #out1 is the final value of our trial solution for the guess Energy=start
    # the element [-1][0] picks the last element in the array.
    out2=odeint(shoot,np.array([1.0,potential(-20,width)-stop]),space,args=(width,stop))[-1][0]
    #out2 is the final value of our trial solution for the guess Energy=stop.
    guess=(stop+start)/2.0
    #guess is the midpoint of the interval.
    out3=odeint(shoot,np.array([1.0,potential(-20,width)-guess]),space,args=(width,guess))[-1][0]
    #out3 is the final value of the trial solution for the guess Energy=guess.
    print(out1,out2,out3) #prints out the three values. Check that two of them have different signs.
    count=1
    while count<50:
        if np.sign(out1)==np.sign(out3):
            start=guess
        else:
            stop=guess
        out1=odeint(shoot,np.array([1.0,potential(-20,width)-start]),space,args=(width,start))[-1][0]
        out2=odeint(shoot,np.array([1.0,potential(-20,width)-stop]),space,args=(width,stop))[-1][0]
        guess=(stop+start)/2.0
        out3=odeint(shoot,np.array([1.0,potential(-20,width)-guess]),space,args=(width,guess))[-1][0]
        print(out1,out2,out3)
        count=count+1
    return guess
```


```python
find_harmonic(0,0.2,10)
```

    5.0201383839734714e+17 -3125212648653843.5 -384102462.1876018
    5.0201383839734714e+17 -384102462.1876018 8.179017955151048e+16
    8.179017955151048e+16 -384102462.1876018 2.2735542460311012e+16
    2.2735542460311012e+16 -384102462.1876018 8410199084143735.0
    8410199084143735.0 -384102462.1876018 3609208269503528.0
    3609208269503528.0 -384102462.1876018 1670918860379850.0
    1670918860379850.0 -384102462.1876018 803802150803844.9
    803802150803844.9 -384102462.1876018 394198592114886.0
    394198592114886.0 -384102462.1876018 195199516651800.16
    195199516651800.16 -384102462.1876018 97127887432922.77
    97127887432922.77 -384102462.1876018 48446233851932.05
    48446233851932.05 -384102462.1876018 24193597310250.742
    24193597310250.742 -384102462.1876018 12089281720052.924
    12089281720052.924 -384102462.1876018 6042618784911.379
    6042618784911.379 -384102462.1876018 3020660087803.171
    3020660087803.171 -384102462.1876018 1510024184243.1016
    1510024184243.1016 -384102462.1876018 754791579975.5231
    754791579975.5231 -384102462.1876018 377196504924.82135
    377196504924.82135 -384102462.1876018 188404540745.5495
    188404540745.5495 -384102462.1876018 94009784066.47076
    94009784066.47076 -384102462.1876018 46812747577.309586
    46812747577.309586 -384102462.1876018 23214300172.8624
    23214300172.8624 -384102462.1876018 11415094143.450294
    11415094143.450294 -384102462.1876018 5515495667.972689
    5515495667.972689 -384102462.1876018 2565696967.655882
    2565696967.655882 -384102462.1876018 1090797496.5233316
    1090797496.5233316 -384102462.1876018 353347578.82097346
    353347578.82097346 -384102462.1876018 -15377490.00776902
    353347578.82097346 -15377490.00776902 168985051.58826664
    168985051.58826664 -15377490.00776902 76803753.95833339
    76803753.95833339 -15377490.00776902 30713151.41413983
    30713151.41413983 -15377490.00776902 7667832.204175376
    7667832.204175376 -15377490.00776902 -3854834.932364304
    7667832.204175376 -3854834.932364304 1906466.8749052437
    1906466.8749052437 -3854834.932364304 -974183.6211656011
    1906466.8749052437 -974183.6211656011 466144.4477593707
    466144.4477593707 -974183.6211656011 -254012.49864809486
    466144.4477593707 -254012.49864809486 106057.8941995735
    106057.8941995735 -254012.49864809486 -74009.43835303216
    106057.8941995735 -74009.43835303216 16031.911087466704
    16031.911087466704 -74009.43835303216 -28984.9512280908
    16031.911087466704 -28984.9512280908 -6473.4396709923
    16031.911087466704 -6473.4396709923 4792.010075300529
    4792.010075300529 -6473.4396709923 -859.5076897088413
    4792.010075300529 -859.5076897088413 1960.345044902494
    1960.345044902494 -859.5076897088413 558.2694034789339
    558.2694034789339 -859.5076897088413 -133.45077659686976
    558.2694034789339 -133.45077659686976 208.10349588346784
    208.10349588346784 -133.45077659686976 67.79327882835169
    67.79327882835169 -133.45077659686976 -30.900100527368945





    0.09999999922386975




```python
energy_guess=0.09999999922386982  #this is a higher order state
start=np.array([1.0,potential(-20,10)-energy_guess])

xout=odeint(shoot,start,space,args=(10,energy_guess))
xt=np.transpose(xout)
psi=xt[0]
psiprime=xt[1]
plt.figure()
plt.plot(space,psi,linewidth=3,label='wavefunction')
plt.plot(space,psiprime,linewidth=3,label='derivative')
plt.legend()
plt.show()
```

It can be seen that this looks like the ground state solution to Schrodinger's equation. It is symmetric about $x=0$ and has a single maximum. Note that as before the solution is un-normalised. This is not a major issue since it can be handled by just dividing the solution by an appropiate factor.

If our method works then we should be able to find the next solution which will be anti-symmetric with $\psi(0)=0$ and have two extrema. So lets try a new guess.


```python
find_harmonic(0.1,0.5,10)
```

    -384102462.1876018 296197.98766152654 -6477297.217795363
    -6477297.217795363 296197.98766152654 59861035832003.46
    -6477297.217795363 59861035832003.46 107060199690817.53
    -6477297.217795363 107060199690817.53 93390716526662.66
    -6477297.217795363 93390716526662.66 60663938464424.32
    -6477297.217795363 60663938464424.32 34436673456375.676
    -6477297.217795363 34436673456375.676 18328664980746.496
    -6477297.217795363 18328664980746.496 9452929477653.213
    -6477297.217795363 9452929477653.213 4800011019612.841
    -6477297.217795363 4800011019612.841 2418569519148.3726
    -6477297.217795363 2418569519148.3726 1213945480628.1665
    -6477297.217795363 1213945480628.1665 608141165405.1974
    -6477297.217795363 608141165405.1974 304359518018.6893
    -6477297.217795363 304359518018.6893 152250177820.7762
    -6477297.217795363 152250177820.7762 76140448337.71689
    -6477297.217795363 76140448337.71689 38071227916.57752
    -6477297.217795363 38071227916.57752 19033511306.56135
    -6477297.217795363 19033511306.56135 9513802022.609184
    -6477297.217795363 9513802022.609184 4753733946.002121
    -6477297.217795363 4753733946.002121 2373646222.050865
    -6477297.217795363 2373646222.050865 1183589052.9639757
    -6477297.217795363 1183589052.9639757 588556982.6836874
    -6477297.217795363 588556982.6836874 291040132.4514977
    -6477297.217795363 291040132.4514977 142281486.06697762
    -6477297.217795363 142281486.06697762 67902149.51525141
    -6477297.217795363 67902149.51525141 30712393.4219151
    -6477297.217795363 30712393.4219151 12117560.27387205
    -6477297.217795363 12117560.27387205 2820055.6336696185
    -6477297.217795363 2820055.6336696185 -1828545.9350546708
    -1828545.9350546708 2820055.6336696185 495808.012813897
    -1828545.9350546708 495808.012813897 -666409.1133511595
    -666409.1133511595 495808.012813897 -85279.39036267347
    -85279.39036267347 495808.012813897 205205.43215363577
    -85279.39036267347 205205.43215363577 59980.558559482015
    -85279.39036267347 59980.558559482015 -12664.666696325607
    -12664.666696325607 59980.558559482015 23651.491564492022
    -12664.666696325607 23651.491564492022 5439.9809942797865
    -12664.666696325607 5439.9809942797865 -3613.4336624834004
    -3613.4336624834004 5439.9809942797865 964.3693145147117
    -3613.4336624834004 964.3693145147117 -1297.3708430636818
    -1297.3708430636818 964.3693145147117 -144.17968917734882
    -144.17968917734882 964.3693145147117 326.56984179683894
    -144.17968917734882 326.56984179683894 63.84344499949393
    -144.17968917734882 63.84344499949393 -46.33244102958592
    -46.33244102958592 63.84344499949393 23.552718176176754
    -46.33244102958592 23.552718176176754 28.110174948148835
    -46.33244102958592 28.110174948148835 -27.35358833472384
    -27.35358833472384 28.110174948148835 16.978776594366185
    -27.35358833472384 16.978776594366185 19.124047594584443
    -27.35358833472384 19.124047594584443 -29.911842306060564





    0.3000000010381232




```python
energy_guess=0.3000000010381232  #this is the second eigenvalue
start=np.array([1.0,potential(-20,10)-energy_guess])

xout=odeint(shoot,start,space,args=(10,energy_guess))
xt=np.transpose(xout)
psi=xt[0]
psiprime=xt[1]
plt.figure()
plt.plot(space,psi,linewidth=3,label='wavefunction')
plt.plot(space,psiprime,linewidth=3,label='derivative')
plt.legend()
plt.show()
```

Again our code has found the solution that we expect. We can see that it is quite a bit wider than the fundamental solution. We can carry on and find the next solution.


```python
find_harmonic(0.4,0.6,10)
```

    59861035832003.46 -1961884369626.1562 296197.98766152654
    296197.98766152654 -1961884369626.1562 -3136916629505.5293
    296197.98766152654 -3136916629505.5293 -2575830852326.186
    296197.98766152654 -2575830852326.186 -1621298470266.7778
    296197.98766152654 -1621298470266.7778 -905660254960.2147
    296197.98766152654 -905660254960.2147 -478126888035.0189
    296197.98766152654 -478126888035.0189 -245585108719.60425
    296197.98766152654 -245585108719.60425 -124448038960.54208
    296197.98766152654 -124448038960.54208 -62640941945.26467
    296197.98766152654 -62640941945.26467 -31424879707.666275
    296197.98766152654 -31424879707.666275 -15738406402.658878
    296197.98766152654 -15738406402.658878 -7875625969.774823
    296197.98766152654 -7875625969.774823 -3939177748.3140407
    296197.98766152654 -3939177748.3140407 -1969870192.504501
    296197.98766152654 -1969870192.504501 -984852819.5584898
    296197.98766152654 -984852819.5584898 -492174324.9449932
    296197.98766152654 -492174324.9449932 -245863802.8279699
    296197.98766152654 -245863802.8279699 -122867211.68592966
    296197.98766152654 -122867211.68592966 -61285939.857209556
    296197.98766152654 -61285939.857209556 -30494985.608838335
    296197.98766152654 -30494985.608838335 -15099416.593094796
    296197.98766152654 -15099416.593094796 -7401623.256598429
    296197.98766152654 -7401623.256598429 -3552702.1473352667
    296197.98766152654 -3552702.1473352667 -1628244.3839049963
    296197.98766152654 -1628244.3839049963 -665985.6262924892
    296197.98766152654 -665985.6262924892 -184901.59478308292
    296197.98766152654 -184901.59478308292 55622.20711500455
    55622.20711500455 -184901.59478308292 -64612.311253855885
    55622.20711500455 -64612.311253855885 -4501.085109296145
    55622.20711500455 -4501.085109296145 25549.490807977607
    25549.490807977607 -4501.085109296145 10561.39042236284
    10561.39042236284 -4501.085109296145 3013.1869732861933
    3013.1869732861933 -4501.085109296145 -762.1380334701074
    3013.1869732861933 -762.1380334701074 1143.031452632145
    1143.031452632145 -762.1380334701074 172.76231706524854
    172.76231706524854 -762.1380334701074 -268.9307216354026
    172.76231706524854 -268.9307216354026 -58.47879851418564
    172.76231706524854 -58.47879851418564 74.68102823835527
    74.68102823835527 -58.47879851418564 4.939478546828507
    4.939478546828507 -58.47879851418564 -32.817590286098074
    4.939478546828507 -32.817590286098074 -6.745527095575517
    4.939478546828507 -6.745527095575517 45.37220841679595
    45.37220841679595 -6.745527095575517 29.88505730526544
    29.88505730526544 -6.745527095575517 -25.15117253520264
    29.88505730526544 -25.15117253520264 -30.954991368633745
    29.88505730526544 -30.954991368633745 1.3504333634671137
    1.3504333634671137 -30.954991368633745 6.536750480126067
    6.536750480126067 -30.954991368633745 -12.430005399896674
    6.536750480126067 -12.430005399896674 2.721509247961896
    2.721509247961896 -12.430005399896674 19.98266863416973





    0.5000000018347055




```python
energy_guess=0.5000000018347055  #this is a higher order state
start=np.array([1.0,potential(-20,10)-energy_guess])

xout=odeint(shoot,start,space,args=(10,energy_guess))
xt=np.transpose(xout)
psi=xt[0]
psiprime=xt[1]
plt.figure()
plt.plot(space,psi,linewidth=3,label='wavefunction')
plt.plot(space,psiprime,linewidth=3,label='derivative')
plt.legend()
plt.show()
```

We know for the harmonic potential the eigenvalues are evenly spaced. In our case the spacing is about 0.2 so we expect the next solution to be at 0.7.


```python
energy_guess=0.7  #no need to use find now since we know the energy spacing.
start=np.array([1.0,potential(-20,10)-energy_guess])

xout=odeint(shoot,start,space,args=(10,energy_guess))
xt=np.transpose(xout)
psi=xt[0]
psiprime=xt[1]
plt.figure()
plt.plot(space,psi,linewidth=3,label='wavefunction')
plt.plot(space,psiprime,linewidth=3,label='derivative')
plt.legend()
plt.show()
```


```python

```


```python

```

0c76baabbaafa1be065f1df67681be6a02b74400
ipynb
Jupyter Notebook
schrodinger1.ipynb
hjelleyman/lorentz
9b21d9c2ece26604a29157191d6d32aa5d6b8823
MIT
schrodinger1.ipynb
hjelleyman/lorentz
9b21d9c2ece26604a29157191d6d32aa5d6b8823
MIT
2020-03-24T17:52:38.000Z
2021-08-23T20:31:03.000Z
schrodinger1.ipynb
hjelleyman/lorentz
9b21d9c2ece26604a29157191d6d32aa5d6b8823
MIT
Qwen/Qwen-72B
1. YES
2. YES