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

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

自洽求解薛定谔方程的系统势阱图

下面是非自洽方法的程序

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))

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