- #1
- 2,802
- 605
Hi
I'm trying to figure out in what energies the bound states of a particle in a symmetric triangular well happen to be, using MatLab(Actually GNU Octave!)
At first I non-dimensionalized the associated Schrödinger equation and got:
<br /> \frac{d^2 \varphi}{dy^2}+[ \varepsilon-v(\pm 2y-1) ]\varphi=0<br />
Where the plus and minus signs refer to two sides of the well.
Then I calculated the classical turning points (where the total energy equals the potential energy) and considered a multiple of it as a measure of the distance that the wave function should go to zero:
<br /> y=\frac 1 2 (\frac \varepsilon v +1)<br />
In this program, I calculate thes stationary states from zero to a multiple of the mentioned point and store its value at that point. Then plot the end points vs the energy to see in what energies the end point is zero.
The problem is, I'm getting NaN as end points! What is wrong?
Thanks
Ohhh...this is the code and I'm using Numerov's method.
I'm trying to figure out in what energies the bound states of a particle in a symmetric triangular well happen to be, using MatLab(Actually GNU Octave!)
At first I non-dimensionalized the associated Schrödinger equation and got:
<br /> \frac{d^2 \varphi}{dy^2}+[ \varepsilon-v(\pm 2y-1) ]\varphi=0<br />
Where the plus and minus signs refer to two sides of the well.
Then I calculated the classical turning points (where the total energy equals the potential energy) and considered a multiple of it as a measure of the distance that the wave function should go to zero:
<br /> y=\frac 1 2 (\frac \varepsilon v +1)<br />
In this program, I calculate thes stationary states from zero to a multiple of the mentioned point and store its value at that point. Then plot the end points vs the energy to see in what energies the end point is zero.
The problem is, I'm getting NaN as end points! What is wrong?
Thanks
Ohhh...this is the code and I'm using Numerov's method.
PHP:
clear
clc
v=3;
P=-1;
de=.01;
e=[-18:de:-17];
m=length(e);
if(P==1)
phi(1)=1;
dphi=0;
else
phi=0;
dphi=1;
end
for j=1:m
CTP=.5*((abs(e(j))/v)+1);
dy=.001;
y=[0:dy:CTP];
n=length(y);
y(1)=0;
phi(2)=phi(1)+dphi*dy+(1/2)*(-v-e(j))*phi(1)*dy^2+(1/6)*(2*v*phi(1)+(-v-e(j))*dphi)*dy^3+(1/24)*(4*v*dphi+((-v-e(j))^2)*phi(1))*dy+(1/120)*(4*v*(-v-e(j))*phi(1)+((-v-e(j))^2)*dphi)*dy^5;
for i=2:n-1
phi(i+1)=(2*phi(i)-phi(i-1)+(dy^2)*phi(i)*(v*(2*y(i)-1)-e)+(1/12)*(dy^2)*(phi(i-1)*(v*(2*y(i-1)-1)-e))-2*phi(i)*(v*(2*y(i)-1)-e))/(1-(1/12)*(dy^2)*(v*(2*y(i+1)-1)-e));
end
ephi(j)=phi(n)
end
plot(e,ephi)