Research Journal of Applied Sciences

Year: 2010
Volume: 5
Issue: 6
Page No. 383 - 387

Suggestion of Proper Boundary Conditions to Solving Schrodinger Equation for Different Potentials by Runge-Kutta Method

Authors : A. Binesh, A.A. Mowlavi and H. Arabshahi

Abstract: We have solved numerically the Schrodinger wave equation for one-dimension or spherically symmetric potential by writing a code in FORTRAN. The method which we have used in the calculation is based on Runge-Kutta method. The Schrodinger equation has been solved for different V (x) and has been checked for the Hydrogen atom. The solutions are consistent with resent data.

How to cite this article:

A. Binesh, A.A. Mowlavi and H. Arabshahi, 2010. Suggestion of Proper Boundary Conditions to Solving Schrodinger Equation for Different Potentials by Runge-Kutta Method. Research Journal of Applied Sciences, 5: 383-387.


The more general 3-D time-dependent Schrodinger's equation is given by Gasiorowicz (1995):


The right-hand side of the equation is the total energy of system where we have the identification that energy is given by:




In the above equations Ê and H are operators and act on the wave function to give the actual energies. If the potential dos not vary with time one can separate the time-dependent Schrodinger equation to two equations, time dependent and time independent equations; therefore the Schrodinger equation is simplified to:


The solution to the Schrodinger equation in three dimensions is quite complicated in general. Fortunately, nature lends us a hand, since most physical systems are rotationally invariant. It means potential depends on the distance between particles but not their direction. In that case one can separate variables by using suitable coordinates. In case the total angler momentum is conserved, the partial Schrodinger equation can be solved by the separation of variables.

One can substitute separated solution into the time-independent Schrodinger equation and obtain the radial equation which depends on one variable. In recent years there has been much activity in the area of the numerical solution of the radial or one-dimensional Schrodinger equation (Avdelas and Simos, 2001). The aim of this activity is to present a method to solve one-dimension time independent Schrodinger equation (Eq. 12) or the radial Schrodinger equation (Eq. 24) for several potentials and attain full eigenfunctions and eigenvalues. We have written a code which has been written in FORTRAN. The numerical method which we have used is the Runge-Kutta method. The method which has been used is suitable when the partial time dependent Schrodinger equation is separated to ordinary time-independent Schrodinger equations.


Many differential equations cannot be solved analytically, in which case we have to satisfy ourselves with an approximation to the solution. The method which we have used in the calculation is the fourth-order Runge-Kutta method. In numerical analysis, the Runge-Kutta methods are an important family of implicit and explicit iterative methods for the approximation of solutions of ordinary differential equations. They are very accurate, powerful and well-behaved for a wide range of problems.

To explain the Runge-Kutta manner, we start with the single variable differential equation y' = f (t; y) with initial condition y (t0) = y0. Suppose that y0 is the value of the variable at point t0. The Runge-Kutta formula takes y0 and t0 and calculates an approximation for y0+1 at a brief point later, t0 + h. It uses a weighted average of approximated values of f (t; y) at several points within the interval (y0; y0+h). The formula is given by:







Thus, the next value (yn+1) is determined by the present value (yn) plus the product of the size of the interval (h) and an estimated slope. The slope is a weighted average of slopes:

k1 is the slope at the beginning of the interval

k2 is the slope at the midpoint of the interval, using slope

k1 to determine the value of y at the point tn + h/2 using Euler's method (Press et al., 1992)

k3 is again the slope at the midpoint but now using the slope k2 to determine the y-value

k4 is the slope at the end of the interval with its y-value determined using k3

In averaging the four slopes, greater weight is given to the slopes at the mid-point:


The RK4 method is a fourth-order method, meaning that the error per step is on the order of h5, while the total accumulated error has order h4 (Press et al., 1992). Note that the above formulas are valid for both scalar and vector-valued functions (i.e., y can be a vector). One-dimensional Schrodinger equation is a second-order differential equation and any second-order differential equation can be written as a set of two first-order equations:


Each of these first-order differential equations can be solved by the above Runge-Kutta forth orders (Lxaru, 1984). We applied this mathematical method for one-dimensional Schrodinger equation to find the eigenfunctions and eigenvalues for different potential functions:



The Schrodinger equation doses not have the complete algebra solution for an arbitrary potential; however, there is a chance to be solved by numerical methods. In general the Schrodinger equation is a partial differential equation in terms of x, y, z and t but in some cases it is separated to ordinary differential equation. When the potential for a two particle system depends on distance between particles or in the special condition it just depends on one variable ( for example r or x), The Schrodinger equation can be separated to Eq. 12 or 24 which is one-dimensional equation. Before solving the Eq. 12 or 24 for different potentials, we should set the boundary conditions in the numerical calculation according to physics problem. As we know the wave function and the derivative of the wave function must be continuous, finite and the integral of the probability density must be limited:


To overcome these conditions in the calculation we set one of the following conditions in the program. We know when x goes infinity, the wave function and the derivative of the wave function should goes zero. Instead of these conditions we presume:





Here is small number and R is very big number. Applying each one of the above conditions is arbitrary; however, choosing suitable condition in the program is very important to avoid of over flowing in the computer program. In the computer program, we also use the atomic units instead of SI units because Plank constant is a small number for a computer program. In the atomic units length is given in terms of Angstrom Å and energy is shown in terms of electron volt (eV), so we have:


and the fine structure is shown by:


Here, me is the rest mass of electron. We have solved numerically the Schrodinger equation for some different potential. Some of the potentials which we have chosen to test the program and the other is correspond to the potential of the Hydrogen atom.


Example I: As a first example we have solved the Eq. 12 for the following potential:


Here, we set:


Because the potential is symmetry, we have used the following boundary conditions:


The wave function goes constant at x = xmin or xmax, so we choose Ψ (x = xmin or xmax) = 1 and as initial conditions. We start with E1 = 1 eV (ground state) as a arbitrary value and calculate the wave function (by solving the Eq. 12 numerically). If the wave function matches with the boundary conditions, it is a solution; otherwise we change the value of E1 to find a wave function which is consistent with the boundary conditions. Figure 1 shows that E1 = 2:047459 eV is a reliable solution because the wave function becomes zero for the large value of x. Fig. 2 shows the solutions of the Eq. 12 for three different values of E1 (2:047 eV; 2:047459 eV; 2:048 eV ). It shows consistent solution with the boundary conditions occurs when E1 equals with 2:047459.

Fig. 1:

The solutions of the Schrodinger equation for different energy values (E1) of grand state. It shows the suitable solution for the equation occurs when E1 = 2:047459

Fig. 2:

The solutions of the Schrodinger equation for three value of E = 2.048, 2.047459 and 2.047 eV

Fig. 3:

The probability density of the grand, the first and the second excited states

Figure 3 plots the probability density for the grand state and the first and the second excited state.

Fig. 4:

The wave functions of the grand and the first excited states

Fig. 5:

The energy spectrum and the difference between the energy levels

Example II: In example II, we have solved the Schrodinger equation with a quite difficult potential:


The potential is not symmetry; therefore we have chosen the following boundary conditions:


Here we have used R = 20 Å. As we have done before, we start with an arbitrary point for En (energy eigenvalue); find the solution of the Schrodinger equation (the wave function). Then, if the solution is consistent with boundary conditions, it is acceptable; otherwise we exam with different value of En. We continue procedure to find the reasonable solution. In Fig. 4, the wave functions of the ground and the first excited state are presented. The energy eigenvalues for grand state and excited states are shown in Fig. 5.

Hydrogen-like atoms: In a hydrogenic atom or ion with nuclear charge +Ze there is the Coulomb attraction between electron and nucleus. This has spherical symmetry potential only depends on r. This is known as a central potential. For a central potential the radial Schrodinger equation is given by:


Here, is the reduced mass:

l = Angular momentum quantum number
R (r) =

The radial wave function. For Hydrogen or electron-positron atom, the potential is given by:


Therefore, the radial Schrodinger equation is given by:


To solve the radial Schrodinger equation with Hydrogen (or electron-positron) atom potential, the following boundary conditions have been chosen:


Then, an arbitrary value for E has been chosen and the equation has been solved numerically to obtain the wave function. If the obtained wave function is consisted with the boundary conditions, it is solution; otherwise we solve the radial Schrodinger equation with different value for E. We continue this way to get correct solution. By this method, we have calculated the energy of the grand state of Hydrogen and electron- positron atom -13:6056 and -6:803 eV, respectively.

Moreover, Fig. 6 shows the wave function of the grand state for Hydrogen and electron-positron atom. In Fig. 7 we have presented the probability density of the grand state.

Fig. 6:

The radial wave functions of hydrogen and positron-electron atoms

Fig. 7:

The probability density of the grand state

In the above examples we have solved one-dimension Schrodinger equation numerically by using fourth-order Runge-Kutta method.


In this study, we have investigated how fourth-order Runge-Kutta method can be applied in the solution of the Schrodinger equation. Suggestion of proper boundary conditions to solving Schrodinger equation for different potentials is very important in numerical solutions. Using a Runge-Kutta based algorithm we have been able to regularize the problems which rise in solving Schrodinger equation such as the singularity and the infinity of the wave function. For different potentials we have considered suitable boundary conditions. For a problem with a potential which behaves as a Coulomb potential both around the origin and in the asymptotic range we considered a more accurate treatment of the numerical boundaries.

Design and power by Medwell Web Development Team. © Medwell Publishing 2023 All Rights Reserved