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

1. Jan 8, 2014

### evamaster

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:

Code (Text):

#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;
}

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.

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

2. Jan 9, 2014

### maajdl

You should give more explanation.
Without it would take too much time to guess what you did based only on your code.
I was surprised by this:

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

I would have expected the wave function to be a complex number.

3. Jan 9, 2014

### Staff: Mentor

Eigenstates of a real potential can be chosen to be purely real.

4. Jan 9, 2014

### evamaster

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.

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.

Yes, http://i.imgur.com/vB9B5.gif

5. Jan 10, 2014

### Staff: Mentor

It seems like the absolute value of the derivative is the same at thepoint where it is discontinuous. Check your PAR factors to see if you don't have too many/too few.

And if I may comment: try to avoid global variables. And why C++? You don't seem to be using any of its functionality.

6. Jan 11, 2014

### evamaster

Why you recomend to avoid global variables? And I use C++ because i have to.