Applications of the Numerov method to simple quantum systems using Python

Numerov's numerical method is developed in a didactic way by using Python in its {\it Jupyter Notebook} version 6.0.3 for three different quantum physical systems: the hydrogen atom, a molecule governed by the Morse potential and for a quantum dot. After a brief introduction to the Numerov method, the complete code to calculate the eigenfunctions and eigenvalues of the hydrogen atom is presented. The necessary code changes to calculate the other two examples are also provided in the sequel.


Introduction
The vast majority of numerical methods, such those of Newton, Euler, Lagrange, Gauss, Fourier, Jacobi, Runge-Kutta and so many others [1], were introduced in the context of applications in physics, astronomy or in other areas of a technical nature, such as aerodynamics. Since then, numerical analysis was not being recognized as a mathematical discipline and this situation persisted during the first four decades of 20th century. Even today, although some numerical methods are taught in physics courses, within the disciplines of mathematics, little emphasis is given to them in physical applications.
In the las decades, the teaching of computing techniques has become more present and also increasingly essential in the development of students from all areas, and therefore, it would not be different for physics teaching. Having said that, it is of paramount importance that we always produce new teaching materials for new technologies such as the Python programming language, which, despite of being relatively new, has already dominated the market to become one of the most important languages today.
We will disclose here a powerful numerical calculation method originally developed by Boris Vasil'evich Numerov [2,3], see also [4,5,6,7], applying it to the time-independent Schrödinger equation describing physical systems like the hydrogen atom, a diatomic molecule governed by the Morse potential and one model for the quantum dot atom [8,9,10]. These three examples will be solved and the parameters needed for each solution using Numerov's method will be shared in their respective sections.
In summary, this work aims to provide the complete code developed in Python with the Jupyter Notebook for the Numerov's numerical method. However, it is important to emphasize that we do not aim to teach Python to the reader, who must have a basic knowledge of programming to be able to keep up the examples.

Numerov's Method
Numerov's initial motivation was to be able to calculate corrections to the trajectory of comet Halley. Therefore, Numerov's method was initially developed to determine solutions to eigenvalue problems associated with ordinary differential equations of second order of celestial mechanics, which did not contain terms involving the first derivative of a function unknown y(x), that is, equations of the form d 2 y dx 2 = f (y, x). (1) Every differential equation equal to equation (1) can be replaced by the following system of first order equations Traditional methods for numerically solving this system of equations, such as those of Euler or Runge-Kutta, consider that the values of y(x) and of dy/dx are known at a given point in the domain [a, b] of system validity, i.e., are suitable for the so-called seed problems.
In non-relativistic quantum mechanics, more specifically in bound state problems involving a particle of mass m confined in a well of potential V (x), in a given interval a < x < b, the allowed energies (E) and the corresponding wave functions ψ(x) that describe these steady states satisfy the Schrödinger's eigenvalue equation In these cases, as the value of the first derivative of the wave function is not known, the Euler and Runge-Kutta cannot be employed. Nonetheless, it is possible to establish continuity conditions for the values of ψ and dψ/dx at two or more points of the domain of the wave function, which characterizes the so-called boundary value problems.
In addition to making the transformation of a second order differential equation in a first order system, the Numerov method allows the simultaneous determination of the energy spectrum of the particle and of the eigenfunctions associated with each energy value.
Like any iterative numerical method, the solution of equation (2) is constructed by successive integrations. In Numerov's method, initially, the solution is considered to be known at two subsequent points of the interval [a, b], for example, at ψ(x − δ) and ψ(x), where δ is an arbitrarily small quantity, called the integration step. Next, we try to establish an algorithm to determine the solution at the next point, ψ(x + δ).
The starting point for establishing this algorithm is the expansion of ψ(x ± δ) in Taylor series, up to fourth-order derivatives, that is, Adding the terms ψ(x + δ) and ψ(x − δ), only the derivatives of even order survive and, therefore, a relationship between the values of a function in three is reached. points and its second derivative, given by Writing the unidimensional Schrödinger equation, equation (2), in a more convenient form and using equation (4) to replace the terms that contain second order derivatives, one obtains Regrouping the therm we obtain the Numerov difference formula for the problem of a particle under action of a one-dimensional potential In fact, it should be noted that the algorithm can be applied to any ordinary linear differential equation and second-order homogeneous that does not contain terms of first derivative.
Since the problem of interest is an eigenvalue problem, the numerical integration technique of the one-dimensional Schrödinger equation for a particle in a well depends on attaching arbitrary values conveniently to eigenvalues and to the respective (possible) eigenfunctions in 2 points of the domain of the problem. But how to do it? Regarding the choice of the initial value for the energy (first eigenvalue), just remember that, according to Heisenberg uncertainty relation, the energy E of a particle in a well of potential V (x) must be greater than the minimum value of the well. Thus, it is considered, initially, that E initial = V min +∆E, with ∆E > 0 The choice of a energy value, determines two turning points, x and x r , where the energy value is equal to potential energy value, whose motion obeys the classical Newtonian mechanics. That is, from the point of view of classical mechanics, the movement of the particle is restricted only to the region [x , x r ], in which the energy is greater than or equal to the potential energy. The regions x < x and x > x r are called classically prohibited regions, and are indicated in Fig. 1 As the Schrödinger equation admits solutions for these classically prohibited regions, for each energy value, initially, values are assigned to a possible eigenfunction at two points of the classically prohibited regions, in which the function practically cancels itself. In general, these are the boundary points a and b > a of the function's integration domain.
However, the implementation of Numerov's method to solve the problem still requires an iteration scheme that uses the Numerov formula in two steps: from a, or to the left of a from the classic turning points, hereinafter called match point (x match ), and from b, or to the right of the match point. Thus, arbitrarily taking a initial value for the energy, and two successive arbitrary values for the solution, starting from the lower extremes and upper part of the integration interval [a, b], one can implement the method's iteration scheme in the two senses, such as: • Solution to the left of match point (x < x match ) Being E initial = V min + ∆E (∆E/|V min | 1) an arbitrary value for the energy of the particle. Also arbitrating values for the function of wave, in 2 successive points, from a, ψ (a) = 0, ψ (a + δ) = δ , (δ 1) and using the formula of differences, equation (7), the solution on the left is built sequentially until match point (x match ), in what ψ (x match ) = ψ match .
• Solution to the right of match point (x > x match ) From a similar way, for the same values E initial , arbitrating the solution to the right, from b, is constituted sequentially until the points x match and x = x match − δ, as To guarantee the boundary condition of the solution, we redefine the solution to the left according to equation (8) given below, and the boundary condition of the first derivatives, according to equation (9).
The procedure is repeated step by step, in the two ways, a b. Starting from a, using the Numerov recurrence formula associated with an equation, if we build the solution ψ until the classic rewind point, nearest of b, where E = V (x match ), called the match point. Then, from b, the analog is made, building a solution ψ r to the match point. In principle, the possible solutions ψ and ψ r will not necessarily be equal in this x match stitch. To ensure the continuity of the solution redefines itself ψ like Finally, it is verified how close are the values of the respective first derivatives of ψ r and the new function ψ It is staggered, at match point. To test the boundary condition of the derivatives first, taking into account Taylor's series for ψ(x + δ and ψ(x − δ), up to the first order, you can write in what, ψ match±1 = ψ(x match±δ ).
If the difference between these values is less than the values of a predefined error, the process is interrupted, confirming the searched eigenvalue and the respective eigenfunction as being If the continuity condition of the derivatives is not satisfied, the value of the energy is increased and we restarted a search for a new value which is really an eigenvalue of the problem, and its respective eigenfunction. The process can be repeated until the desired number of eigenvalues and eigenfunctions of the problem.
Because it is based on Taylor's serial expansion to fourth order, the error in Numerov's method is much smaller than the errors that come out from the expansion-based methods in lower order, like that of Runge-Kutta.

Hidrogen Atom
Numerical solutions of hydrogen atom was previously obtained in [1] using Numerov's method. The program was written in C++ for the ROOT cint compiler.
Although it was originally developed for second order linear and homogeneous ordinary differential equations that do not contain terms of the first derivative, the Numerov's method can be generalized to cover the presence of terms that contain the first derivative in the differential equation, so that eigenvalue problems can also be considered.
In fact, in the case of linear equations, every equation second order differential of type can be written in its normal form Schrödinger's radial equation for a particle of mass m under the action of a Coulombian electric field, like the electron in the hydrogen atom, can be written as Making the substitution r = xa B with a B = 2 /(me 2 ) being the Bohr radius, equation (10) can be rewritten, for a new function y(x) = R(r) as x are, respectively, the energy and the so-called effective potential in atomic units. So, in possession of equation (11), we can start building our program code.
First of all, we must import the functions available in the Pylab module that bulk imports matplotlib.pyplot (for plotting) and NumPy (for Mathematics and working with arrays) in a single name space. We then declare who our effective potential V (x) is, and during the construction of this example we will use = 0.
if j<int(N) and abs(x[j])>xmax: The equation, in this case, that is intended to be solved by Numerov's method presents a term involving the first derivative, and can be expressed by where From the Taylor expansions, equation (3), we can rewrite equation (12) as In a similar way to the previous case, according to equation (4), you can write the term on the right side of the equation (13) which contains derivatives of order 2, such as Replacing first order derivatives with approximations Taking into account that the left side of the equation (13) is equal to Regrouping the terms, and making one obtains the Numerov difference equation for the problem, suitable for the propagation of the solution from of the limits of the integration interval From this formula, a procedure analogous to the previous case can be implemented for the construction of solutions of the radial Schrödinger equation in the interval (0, ∞). Now, in order to introduce the Numerov diference formula (15), we first need to insert equation (7) in our code, for that, let's break it down into different pieces p 0 , p 1 and p 2 , with h = δ, q 0 = s(x). Thus, equation (15) is now called print(" x0 = ", x0," y0 = ", y0," V = ",V(x0)) print(" x1 = ", x0+h," y1 = ",y1," V = ",V(x0+h), " y2 = ",y2) return y2 And, for the case x − δ, we repeat the previous step, changing the appropriate sign.
xx=list(range(dim)); yy=list(range(dim)) ww=list(range(dim)); yl=list(range(dim)) yr=list(range(dim)); ee=list(range(nmax)) ff=list(range(nmax)); ff2=list(range(nmax)) yy1=list(range(dim)); yy2=list(range(dim)) yy3=list(range(dim)) colors =['b','r','g','m','c']; nk=list(range(nmax)) E_old = Ein; E = Ein + dE In the next steps, the program will determine the interaction for the eigenvalue candidates and determine their solutions as well as printing the parameters found on the screen.  So far, the only thing we need to change in our program is the equation that defines our effective potential V (x). The rest of the code will be the same for any equation that is in the same form as equation (9). From now on, we must change the code whenever we are looking for solutions with a different potential.
The parameter M indicates the eigenvalue that we are determining, in the next step, in order to avoid that the program needs to sweep the entire potential well in search of solutions, we can give increments between one eigenvalue and another in the form of multiples of the dE parameter, which will be duly defined in a next step.
Usually, this process involves trial and error, where we adjust the dE multiplier for each case, until we find the desired eigenvalue. However, as we are developing this first example for the hydrogen atom, which already has its energy eigenvalues well defined, this task becomes much easier. Before we find the ground state (M = 0) we must use the step 24 × dE. Thus, we find for ground state energy the value of ε = −0.99453, which must be compared with the known value of the ground state energy for the hydrogen atom ε = −1. And the next multipliers so that we can find at least the first three solutions witch are: We finally reach the part of the program where we must introduce the initial values to proceed with the solution. Whenever we start to develop a new code, these must be the first values to be changed, right after choosing the effective potential V (x).
The parameters a and b correspond to the boundary points mentioned in the second section. At first, our wave function should have its initial value (a) equal to 0, however, to avoid divisions by 0 throughout the program we should start from a relatively small value (a = 0.001). In the piece of code below h represents the δ increment, kmax is the maximum number of interactions for each eigenvalue, nmax is the number of eigenvalues that we are trying to find, Ein is the minimum energy of de effective potential that we are trying to solve, and V max is the maximum energy of the same effective potential.

Morse Potential
The Morse potential is a common model for the interatomic interaction of a diatomic molecule [11,12]. In this section, in order to learn how we can use the Phyton code in other problems, we will see what we must change in the code that was made available in the previous section so that the program will be able to solve equation (2) for the quantum number = 0 and the Morse potential with arbitrary parameters given by: In this example, we will calculate the first two bounded states, for that, we must adjust the dE multiplier for each case as: And, the most important part, which is the adjustment of the initial data of our problem. Analyzing the effective potential we can verify that the value os parameter a must be negative, and its minimum value is zero, so the code must be set as follows: From now on, the program is able to find the energies of the first two states, which are respectively 7.1380 and 15.0380. As well as already determined the wave functions also for the first two states.
Then, all that remains is to get all the aesthetic part of the code right, adjusting the limits of graphics, subtitles and other factors. Thus we first get the effective potential, as shown in Figure 6 and finally arrive at the normalized wave functions shown in Figure 7.  In Table 2, we can compare the eigenvalues found by the Numerov's method with the analytical values given by E n = 16 n +

Quantum Dot
The development of technology based on quantum dots is quite recent, but it is already showing signs that it is the next great technology, when we talk about optics. In a simple model for a quantum dot composed of two electrons, they can be described with a external harmonic oscillator potential of frequency Ω = 2ω. Following the steps of reference [10], we have the effective potential to be introduced into equation (2) for the quantum number = 0 is given by: introducing this potential into our code, let's now calculate the first five solutions that the program is capable of finding. We can observe that in our code used so far, we only have up to three wave functions, in this example we will show how we included two more solutions in the code. The procedure is very simple, just add the yy4 and yy5 functions along the code and all the other parts that are related to them, for example, xx=list(range(dim)); yy=list(range(dim)) ww=list(range(dim)); yl=list(range(dim)) yr=list(range(dim)); ee=list(range(nmax)) ff=list(range(nmax)); ff2=list(range(nmax)) yy1=list(range(dim)); yy2=list(range(dim)) yy3=list(range(dim)); yy4=list(range(dim)) yy5=list(range(dim)) colors =['b','r','g','m','c']; nk=list(range(nmax)) E_old = Ein; E = Ein + dE Now, as in the previous example, let's include the dE multipliers, remembering to include the new wave functions. And finally, we must include the initial conditions for our problem, as shown below.
Thus, adjusting the rest of the code, we obtain the first five eigenvalues. The effective potential and the eigenfunctions can be seen in Figures 8 and 9.
In Table 3, we can compare the eigenvalues found by the Numerov's method with the analytical values for the quantum dot given by η n = 2(n + + 1)ω.

Conclusion
Analyzing the results arranged in Tables 1, 2 and 3 we can conclude that the method used here is able to reproduce the analytical results within small errors. Therefore, it is evident that the numerical method of Numerov is a powerful tool, easy to use, which can help in the development of not only new knowledge in the area of programming, but also the solution of Schrödinger equations outside the usual results found in the examples of modern physics books. We hope that this short introduction to the code, which can be found in full at https://1drv.ms/u/s!Ai_Lqkgh1kiskp5_cfOtwCfqX-LpTw?e=srdmzS, will open doors for students to create their own versions, increasingly improving the versatility of this tool.