# Implementation of the Numerov Method for the 1D square well

I want to solve the Schrodinger via the Numerov Method but I had some troubles. I'm programing in C++, so here is my code:

#include<cstdlib>
#include<iostream>
#include<cmath>

using namespace std;

double x_min=-4.0 , x_max=4.0;
int N=2000;
double r=(x_max-x_min)/(1.0*N);
double d=2.0;
double p=0.4829;    // 2m/(hbar^2)
double Vo=20.0; // Altura del pozo

double x_m=0.1;                                     //Matching point
int i_x_m=(x_m-x_min)/r;

double Control=-123456789;

double SlopeLeft,SlopeRight;

double PAR;

double K2(double x, double E);
double NumerovL(int i, double k21, double k22, double k23, double Y[]);
double NumerovR(int i, double k21, double k22, double k23, double Y[]);
double FuncLeft(double E, double Y[]);
double FuncRight(double E, double Y[]);
void PrintFunc(double Y[]);
void Normalizar(double Y[]);
double f(double E, double Y[]);
double Biseccion(double a, double b, double Y[]);

//=========================MAIN===============================

int main(int argc, char **argv)
{

double Y[N+1];       // Función de Onda

double paso=0.02;   // Escala   en la que se varia la energía
double Eo=0;

for(double E=0 ; E<=Vo ; E+=paso) //Cálculo de las funciones IMPARES
{
PAR=-1;
Eo=Biseccion(E,E+paso,Y);

if(Eo != Control && SlopeRight*SlopeLeft<0.)
{
Y[i_x_m]=FuncRight(Eo,Y);
Y[i_x_m]=FuncLeft(Eo,Y);

Normalizar(Y);

PrintFunc(Y);
}

}

for(double E=0 ; E<=Vo ; E+=paso)   //Cálculo de las funciones PARES
{
PAR=1;
Eo=Biseccion(E,E+paso,Y);

if(Eo != Control && SlopeRight*SlopeLeft>0.)
{
Y[i_x_m]=FuncRight(Eo,Y);
Y[i_x_m]=FuncLeft(Eo,Y);

Normalizar(Y);

PrintFunc(Y);
}
}

return 0;
}

//=========================FUNCIONES===============================

double K2(double x, double E)
{
double k2;

if(fabs(x)<=d)
{
k2=p*E;
return k2;
}
else
{
k2=p*(E-Vo);
return k2;
}
}

double NumerovL(int i, double k21, double k22, double k23, double Y[])
{ // Para la función de Onda Izquierda
double A1,B1,C1,N;
A1=2.0*(1.0-(5.0/12.0)*r*r*k21)*Y[i-1];
B1=(1.0+(1.0/12.0)*r*r*k22)*Y[i-2];
C1=1.0+(1.0/12.0)*r*r*k23;
N=(A1-B1)/(C1);
return N;
}

double NumerovR(int i, double k21, double k22, double k23, double Y[])
{ // Para la función de Onda Derecha
double A1,B1,C1,N;
A1=2.0*(1.0-(5.0/12.0)*r*r*k21)*Y[i+1];
B1=(1.0+(1.0/12.0)*r*r*k22)*Y[i+2];
C1=1.0+(1.0/12.0)*r*r*k23;
N=PAR*(A1-B1)/(C1);
return N;
}

double FuncLeft(double E, double Y[])
{
double k21,k22,k23,Yleft,b;

b=sqrt(p*(Vo-E));

Y[0]=exp(b*x_min);
Y[1]=exp(b*(x_min+r));

for(int i=2 ; i<i_x_m ; i++) // Se calcula la función de Onda Izquierda
{
k21=K2(x_min+(i-1)*r,E);
k22=K2(x_min+(i-2)*r,E);
k23=K2(x_min+i*r,E);

Y[i]=NumerovL(i,k21,k22,k23,Y);

if(i==i_x_m-1) //Función de Onda Izquierda en el Matching point
{
k21=K2(x_min+(i)*r,E);
k22=K2(x_min+(i-1)*r,E);
k23=K2(x_min+(i+1)*r,E);

Yleft=NumerovL(i+1,k21,k22,k23,Y);
}
}

SlopeLeft=(Yleft-Y[i_x_m-1])/r;

return Yleft;
}

double FuncRight(double E, double Y[])
{
double k21,k22,k23,Yright,b;

b=sqrt(p*(Vo-E));

Y[N]=PAR*exp(-b*(x_min+N*r));
Y[N-1]=PAR*exp(-b*(x_min+(N-1)*r));

for(int i=N-2 ; i>i_x_m; i--) // Se calcula la función de Onda Derecha
{
k21=K2(x_min+(i+1)*r,E);
k22=K2(x_min+(i+2)*r,E);
k23=K2(x_min+i*r,E);

Y[i]=PAR*NumerovR(i,k21,k22,k23,Y);

if(i==i_x_m+1) //Función de Onda Derecha en el Matching point
{
k21=K2(x_min+(i)*r,E);
k22=K2(x_min+(i+1)*r,E);
k23=K2(x_min+(i-1)*r,E);

Yright=NumerovR(i-1,k21,k22,k23,Y);
}
}

SlopeRight=PAR*(Y[i_x_m+1]-Yright)/r;

return Yright;
}

void PrintFunc(double Y[])
{
for(int i=0 ; i<=N+1 ; i++)
{
cout << x_min+i*r << "\t" << Y[i] << endl;
}
}

void Normalizar(double Y[])
{
double S=0;

for(int i=0 ; i<=N+1 ; i++)
{
S += Y[i]*Y[i]*r;
}

S=sqrt(S);

for (int i=0 ; i<=N+1 ; i++)
{
Y[i]=Y[i]/S;
}

}

double f(double E, double Y[])
{
double F;

F=FuncLeft(E,Y)-PAR*FuncRight(E,Y);

return F;
}

double Biseccion(double a, double b, double Y[])
{

double Tol=0.00001; //Tolerancia para encontrar la raiz

double RET=-123456789;

if(f(a,Y)*f(b,Y)<0)
{
while(fabs(a-b)>Tol)
{
double x_m,fa,fm;

fa=f(a,Y);
x_m=(a+b)/2.0;
fm=f(x_m,Y);
//fb=f(b);

if(fa*fm<0)
{
b=x_m;
//RET=b;
}
else
{
a=x_m;
//RET=a;
}
}
RET=a;
}

return RET;
}


Basically the code takes all the energies, i.e. $0<E<Vo$ and the fuction "Biseccion" applies the Bisection algorithm between an energy $E$ and $E+step$. So the fuction finds the eigen-energy for wich the left and right (from the Numerov Method) wave functions matches.

The code compiles perfectly but the problem arises when I want to plot the odd solutions. I obtain two satisfactory solutions but another two that it's function is continous but not it's derivative. Here is an example of the plot that I obtain:

As you can see, there are two graphs that are not a satisfactory solution to the problem.

I would be very thankful if somebody can help me with this problem.

• It would help to see plots of the solutions in question. Jan 9 '14 at 1:56
• If you wanna check out, I've updated my question with a graph that I obtain with the data from my program. Jan 9 '14 at 3:20
• This doesn't help your specific problem but if you just want to solve the Schrodinger equation the easiest way I have found is to discretise the differential operator (e.g. using finite difference) and simply solve the resulting matrix (a tridiagonal system if you use centred differences) for the eigenvalues and vectors. This real space approach has some draw backs because it may return spurious wavefunctions but they are normally easily filtered because they are clearly non-physical. Can post more info if useful. Jan 9 '14 at 14:03
• Thats a good idea. But I have to solve the problem using ONLY the Numerov Method link . Correct me if I'm wrong, but I think that your solution doesn't apply this method. Because if your solution contains this method, your help would be very useful.Thanks. Jan 9 '14 at 16:16
• Using the method of manufactured solutions in concert with a unit testing framework will help you diagnose if the error is in your implementation of the Numerov method. Jan 9 '14 at 19:43

## 1 Answer

I think @Ondřej-Čertík already pointed it out: I think you are obtaining the right solutions, but notice that the matching condition between the left and right solutions holds even if you multiply one or both solutions for a constant, so you are free to scale one of your functions as to exactly match the other in the matching point, this is equivalent to change the next-to-boundary initial values. Then you can normalize if you finally need it. In any case I think this is enough as to obtain the right eigenvalues.