N=200 %%要素を分ける
V=0.3 %%potential
a1=100 %%potential range
a2=45 %%potential range
a3=65 %%potential range
a=6.2415*10^18*(6.6*10^(-34)/(2*pi))^2/(2*0.067*9.1*10^(-31))
l=(a1/N)*10^(-9)
for i=1:N+1
if i <= a2*(N/a1)
A(i)=1;
else
if i <= a3*(N/a1)
A(i)=0;
else
A(i)=1;
end
end
end
B(N+1,N+1)=0;
for m=1:1:N
BB(N+1,N+1)=0;
BB(m,m)=a * (1/l) + (l/3) * V * A(m);
BB(m,m+1)=a * (-1/l) + (l/6) * V * A(m);
BB(m+1,m)=a * (-1/l) + (l/6) * V * A(m);
BB(m+1,m+1)=a * (1/l) + (l/3) * V * A(m);
B=B+BB;
BB=0;
end
B1=B;
B(N+1,:)=B1(2,:);
B(2,:)=B1(N+1,:);
B2=B;
B(:,N+1)=B2(:,2);
B(:,2)=B2(:,N+1);
B=B(3:N+1,3:N+1);
for m=1:N-2
BB=B;
B(N-1,:)=BB(m,:);
B(m,:)=BB(N-1,:);
BBB=B;
B(:,N-1)=BBB(:,m);
B(:,m)=BBB(:,N-1);
end
C(N+1,N+1)=0;
for m=1:N
CC(N+1,N+1)=0;
CC(m,m)=l/3;
CC(m,m+1)=l/6;
CC(m+1,m)=l/6;
CC(m+1,m+1)=l/3;
C=C+CC;
CC=0;
end
C1=C;
C(N+1,:)=C1(2,:);
C(2,:)=C1(N+1,:);
C2=C;
C(:,N+1)=C2(:,2);
C(:,2)=C2(:,N+1);
C=C(3:N+1,3:N+1);
for m=1:N-2
CC=C;
C(N-1,:)=CC(m,:);
C(m,:)=CC(N-1,:);
CCC=C;
C(:,N-1)=CCC(:,m);
C(:,m)=CCC(:,N-1);
end
D=C\B;
[v,e]=eig(D);
e=diag(e);
x=[a1/N:a1/N:a1-a1/N];
xV=[0:a1/N:a1];
clg;
hold on;
for m=1:2
plot(x,v(:,m));
endfor
A=V*A;
plot(xV,A);