Trạng thái dừng trong hố thế vuông

function standing_wave_square_potential_well_new
% Created by Tran Hai Cat
% 2018.05.15
clc;
clear variables
close all

global E k U0 x1 x2

%% CONSTANTS
h = 1.054e-34;
m = 9.1095e-31;
q_element = 1.6e-19;
Angstrom = 1e-10;
k = 2*m/h/h;

%% INPUT PARAMETRS
n = 2; % n<=6, index of stationary state

%% DATA PROCESSING
% Energy of stationary state:
E_standing_wave = [-78.097995,-72.369,-62.905,-49.8624,-33.7612,-14.38695];

E = E_standing_wave(n)*q_element;
U0 = -80*q_element;
x1 = -2*Angstrom;
x2 = 2*Angstrom;

Emin = -100*q_element;
Emax = 100*q_element;
xmin = -6*Angstrom;
xmax = 6*Angstrom;
x_inter = 0.5769*Angstrom;

psi = 1e-10;
psi1 = 0;

%% CALCULATION
[x_left,psi_left] = ode45(@syst,linspace(xmin,x_inter,500),[psi psi1]);
[x_right,psi_right] = ode45(@syst,linspace(xmax,x_inter,500),[psi psi1]);
psi_right = sign(psi_left(end,2)*psi_right(end,2))*psi_right;
X = [x_left;flipud(x_right)];
Psi = [psi_left(:,1);flipud(psi_right(:,1))];
Psi_max = max(Psi);
Psi_draw = 0.7*Psi/Psi_max*Emax/q_element;
psi_psi = 0.8*(Psi.^2/Psi_max/Psi_max*Emax/q_element);

omega = abs(E)/h;
T = 2*pi/omega;
dt = T/1000;
t = 0;

%% FIGURE
x_plot = linspace(xmin,xmax,1000);
U_plot = U_square(U0,x1,x2,x_plot);

figure('name','Standing_wave_potential_well',...
'color','black','numbertitle','off');
hold on
plot(x_plot/Angstrom,U_plot/q_element,'linewidth',2,'color','w');
line([xmin/Angstrom,xmax/Angstrom],[E/q_element,E/q_element],...
'linewidth',2,'color','b');
line_psi = plot(X/Angstrom,Psi_draw,'linewidth',1,'color','y');
plot(X/Angstrom,psi_psi,'linewidth',1);
title(sprintf('E = %0.3f [eV]',E/q_element),'color','w');
axis([xmin/Angstrom xmax/Angstrom Emin/q_element Emax/q_element]);
set(gca,'color','k','xcolor','w','ycolor','w')
xlabel('x [Angstrom]');
ylabel('E, U [eV]');

while 1
t = t+dt;
Psi_time = Psi_draw*exp(-1i*omega*t);
set(line_psi,'ydata',real(Psi_time));
pause(0.001);
end

function dy = syst(x,y)
%% Schrodinger Equation
global E k U0 x1 x2
dy = zeros(2,1);
dy(1) = y(2);
dy(2) = -k*(E-U_square(U0,x1,x2,x))*y(1);

function y = U_square(U0,x1,x2,x)
%% Potential Energy
lenx = length(x);
y = zeros(1,lenx);
for i = 1:lenx
if(x(i)<=x1)||(x(i)>=x2)
y(i) = 0;
else
y(i) = U0;
end
end