C ++: Реализация метода Нумерова для одномерной квадратной скважины

Я хочу разгадать Шредингера по методу Нумерова, но у меня были некоторые проблемы. Я программирую на C ++, поэтому вот мой код:

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

В основном код забирает все энергии, т.е. 0<E<Vo и функция «Biseccion» применяет алгоритм деления пополам между энергией E а также E+step, Таким образом, функция находит собственную энергию, для которой левый и правый (из Метод нумерова) волновые функции совпадают.

Код компилируется отлично, но проблема возникает, когда я хочу построить нечетные решения. Я получаю два удовлетворительных решения, но еще два, что его функция непрерывна, но не является ее производной. Вот пример участка, который я получаю.

Как видите, есть два графика, которые не являются удовлетворительным решением проблемы.

Я был бы очень благодарен, если кто-то может помочь мне с этой проблемой.