# One Dimensional Schrodinger's Equation solution using Numerov Method

I have been trying to solve Time Independent Schrodinger's equation in one dimension using Numerov Method as discussed in this excellent lecture notes I found on net.

The Numerov method can solve an equation of the following kind:

$$\frac{{d^2}y}{dx^2}=-g(x) y(x) +s(x)$$

We can compare this with out Time Independent Schrodinger Equation :

$$\frac{{d^2}\psi}{dx^2}=-\frac{2m}{\bar h^2} \big(E-V) \psi$$

Where $g(x)= -\frac{2m}{\bar h^2}(E-V)$

We can discretize the equation using Taylor expansion and write the iterative formula as :

$$\psi_{i+1}= \frac{\psi_i\big[2-\frac{5}{6}\cdot h^2 g_i]-\psi_{i-1}\big[1+\frac{h^2}{12}\cdot g_{i-1}] }{(1+\frac{h^2}{12}\cdot g_i)}$$

Using the known Energy of Harmonic Oscillator I solved the equation and my solutions were working well.

I also used the shooting method to determine Energy value for a particular Eigen State and got a great result.

I was succesfully able to solve the schrodinger equation for Most of Even potential including Morse Potential and The Double Well Potential.

The Problem arrises if I try to take a non symmetric solution. Because I have taken the symmetry of my wave function into account and just solved for $\psi$ from $x=0$ to $x=xmax$, and then took output in a file by multiplying a parity to the $\psi$ based on weather my function is odd or not.

But When I solved the same problem from $x=-xmax$ to $x=xmax$ using same algorithm so that I could devise a general algorithm for Numerov Method. The output started getting bizarre. While the ground state solution is same , the odd state solutions are getting inverted somehow.

This is a part of my program for the first (working) Numerov Algorithm. (ignore the signchange counter I have used it to determine number of nodes for shooting method)

cout<<"Enter mesh number"<<endl;
cin>>n;

cout<<"Enter maximum value of x"<<endl;
cin>>xmax;
dx=xmax/n;

if (nodes % 2==0) //Even Nodes, function is even function , i.e. y(-1) =      y(1) and f(-1) = f (1); Thus using algorithm!
{
y[0]=1.0;
y[1]=( y[0]*( 12.0 - 10.0 *f[0] ) )/( 2*f[1] );
}
else //nodes are odd, i.e there is a node at x=0;
{
y[0]=0;
y[1]=dx;           //arbitrary small value
}

signchange=0;

//outward integration using algorithm and counting number of sign changes//

for(int i=1; i <= n; i++)

{
y[ i+1 ]= ( ( 12.0 - 10.0*f[i] )*y[i] -( y[i-1]*f[i-1] ) ) / f[i+1];

if (y[i] != copysign(y[i],y[i+1]))
++signchange;
}
//For Output :

if (nodes % 2 == 0)
p=1;
else
p=-1;
for (int i = n; i > 0; --i)
file<<setw(1)<<-1*x[i]<<setw(15)<<p*y[i]<<setw(15)<<y[i]*y[i]<<setw(15)<<V[i]<<endl;

//-------x<0-------//
for (int i = 0; i <= n; ++i)
file<<setw(1)<<x[i]<<setw(15)<<y[i]<<setw(15)<<y[i]*y[i]<<setw(15)<<V[i]<<endl;


And this is part of the second program which is giving me inverted output. Where I am iterating across the whole mesh

cout<<"Enter maximum value of x"<<endl;
cin>>xmax;

xmin=-1*xmax;

h=(2*xmax)/n;

y[0]=0;
y[1]=h;

/*calculating the wave function psi or y at all points for Energy e*/
for(int i=1; i <= n; i++)

{
y[ i+1 ]= ( ( 12.0 - 10.0*f[i] )*y[i] -( y[i-1]*f[i-1] ) ) / f[i+1];

if (y[i] != copysign(y[i],y[i+1]))
++signcount;
}


I would love to know why this is happening. It does not make sense to me. But i am also a beginner in Computational Physics. Also I would like to know what should be the best way to solve the equation for Non-Symmetric potentials, for example a Coulumb potential (inverse R).

EDIT: Since I am trying to solve Harmonic Oscillator problem the wave function should be bound on both sides. So i am taking $\psi_0$ as 0. And $\psi_1$ as some random small value which will be taken in to account when the wave function will be normalized later. This is for the case when I am solving the equatioin across the whole mesh, from $-xmax$ to $xmax$ .

But When I am solving the equation for symetric potential(i.e. from 0 to xmax), I am using a property of solution that I have already known, that if i am solving for odd Energy state, the wave funcion $\psi_{x=0}$ will always be zero. I have provided my declaration of initial value in code too. I hope this helps, in clarifying my problem a little.

• Please give a more complete description of your problem by including boundary conditions. You can edit the question to add this information. – cfdlab Sep 7 '17 at 17:07
• What are the lecture notes you are refering to? – nicoguaro Sep 7 '17 at 18:40
• fisica.uniud.it/~giannozz/Corsi/MQ/LectureNotes/mq-cap1.pdf these are the lecture notes I am talking about – pennywise666 Sep 7 '17 at 23:12
• So, what is your problem? The inverted sign is not a problem for eigenfunction calculation. – nicoguaro Sep 8 '17 at 1:19
• It's really not clear here what the problem is supposed to be. Could you perhaps post some kind of an MCVE that shows exactly what answer is being computed and why it is wrong? – Kirill Sep 8 '17 at 4:52