Computer Physics Communications 36 (1985) 113—119 113 North-Holland, Amsterdam A VARIABLE STEP METHOD FOR THE NUMERICAL INTEGRATION OF THE ONE- DIMENSIONAL SCHRODINGER EQUATION AD. RAPTIS Department of Mathematics, National Technical University of Athens, Zografou Campus, Greece and J.R. CASH Department of Mathematics, Imperial College, South Kensington, London S W7 2A Z, England Received 12 November 1984 Most numerical methods which have been proposed for the approximate integration of the one-dimensional SchrOdinger equation use a fixed step length of integration. Such an approach can of course result in gross inefficiency since the small step length which must normally be used in the initial part of the range of integration to obtain the desired accuracy must then be used throughout the integration. In this paper we consider the method of embedding, which is widely used with explicit Runge—Kutta methods for the solution of first order initial value problems, for use with the special formulae used to integrate the Schrodinger equation. By adopting this technique we have available at each step an estimate of the local truncation error and this estimate can be used to automatically control the step length of integration. Also considered is the problem of estimating the global truncation error at the end of the range of integration. The power of the approaches considered is illustrated by means of some numerical examples. 1. Introduction Many problems in application areas such as theoretical chemistry, astrophysics, laser physics, pollution of the atmosphere, nuclear reactions, etc., reduce to a system of second order ordinary differential equations of SchrOdinger form in one independent variable — namely the radius r. The need for improved numerical algorithms for the solution of these equations has been pointed out by many authors (see e.g. ref. [1]). The one-dimensional equation can be considered as being representative of the general problem, which involves coupled systems of equations, and consequently most authors restrict their attention to the study of the scalar equation y”(r) = [v(r) —E] y(r) mf(r)y(r). (1.1) The two boundary conditions associated with (1.1) are y(O) = 0 together with a second condition imposed at large r. The form of the second boundary condition depends on the sign of E, if E> 0 we have a phase shift problem while if E < 0 the problem is of eigenvalue type. The function v(r) appearing in (1.1) is called the potential energy function or simply the potential. In this paper we too restrict our attention to the scalar equation (1.1) but we will be careful to develop techniques which are immediately applicable to vector equations as well. Due to the fact that (1.1) has been studied so extensively, many numerical methods have been proposed for its solution. In ref. [2] at least 20 commonly used methods are cited for the integration of (1.1). One of the most popular of these is Numerov’s method which is the optimal linear two step method [3]. Although Numerov’s method is only of order 4 it has been found in practice to be often superior to higher order linear multistep methods [4]. 0O1O-465~/85/$03.3O © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)