【求助】强者进!!自洽薛定谔--泊松方程解

[ljb0440502112]

自洽求解薛定谔方程的系统势阱图
下面是非自洽方法的程序
clear;
clc;
s0=-250:1:250;
v=0;
j=1;
m=1;
t1=zeros(1,1000000);
t2=zeros(1,1000000);
for E=25:.0001:30
s1=0:1:500;s1(1)=0;s1(2)=1;
for k=2:1:500
if k>=201 & k<=300
v=0;
else v=167;
end
s1(k+1)=(1.9294e-5*(v-E)+2)*s1(k)-s1(k-1);
end
t1(m)=s1(501);t2(m)=E;
m=m+1;
end
for i=2:999999
if abs(t1(i))<abs(t1(i-1))&abs(t1(i))<abs(t1(i+1))
E=t2(i);
s1=0:1:500;s1(1)=0;s1(2)=1;
for k=2:1:500
if k>=201 & k<=300
v=0;
else v=167;
end
s1(k+1)=(1.9294e-5*(v-E)+2)*s1(k)-s1(k-1);
end
h=0;a=0;
for m=1:501
h=h+s1(m)^2;
if h~=0
a=(1/h)^0.5;
end
end
b=a*s1(501);
s1=a*s1;
%s1=s1+.005*E;
subplot(2,2,1)
plot(s0,s1);
e(j)=E;j=j+1;
else continue
end
end
for i=1:501
if i>=201 & i<=300
r1(i)=-(2e16*(s1(i).^2)-2e14)*1e-14;
else
r1(i)=-2e16*(s1(i).^2)*1e-14;
end
end
subplot(2,2,2);
plot(s0,r1);
r2=0:500;
%r2(1)=-0.904*sum(r1(2:501))/13;
for i=1:1:500
r2(i)=.904*(sum(r1(1:i))-sum(r1((i+1):501)))/13;
end
r2(501)=.904*sum(r1(1:501))/13;
subplot(2,2,3);
plot(s0,r2);
%r2=r2-1.808;
r3=0:500;r3(1)=.1*r2(1);
for j=2:501
r3(j)=.1*sum(r2(1:j));
end
subplot(2,2,4);
plot(s0,r3);
%end