摘要：此文介绍了一种使用MATLAB求解一维定态薛定谔方程的方法。利用充分格式进行离散化，得出相应的矩阵方程，用MATLAB求解本征值和本征函数。此方法简单可靠，可以处理各种时间无关的束缚态问题。所用的程序附于文后。
宏观物体遵循的是牛顿运动方程，给定初始条件以及受力条件，我们就可以求出任意时刻粒子的状态。而在原子尺度上，所有粒子都表现出波的行为，粒子的状态用波函数
来描述，描述微观粒子运动的方程是由薛定谔于1925年提出薛定谔方程。微观粒子的运动与其所处的势场相关，当势场不随时间变化时，称为定态，一维情况下的定态薛定谔方程为
其中，
表示粒子所处的势场。由定态波函数可以得出总波函数
1926年， Born提出，波函数模的平方
代表在位置 ，时刻 寻找到电子的概率，这就是波函数的条统计解释。可以由遵循微分方程的波函数表示，薛定谔波方程涉及空间坐标和时间。对于一维情况，在时刻 时，在 和 之间的某处找到电子的概率由下式给出
时间无关薛定谔方程(1)是系统的能量本征方程。该特征值方程通常由一组特定的函数
和相应的一组常数 满足解的条件，被称为是哈密顿算符的本征函数和相应的本征值。当测量处于 状态的系统的能量时，结果将始终为 。对于束缚态情况，必须有
一维定态薛定谔方程是二阶微分方程，但是，能够解析求解的情况屈指可数，如氢原子，谐振子，无限深势阱等。随着计算机技术的发展，我们可以数值上进行求解。本文利用MATLAB软件，使用矩阵方法求解束缚态的本征值问题。
对于原子系统，以nm和eV的能量测量长度更方便。我们可以使用缩放因子
因此我们可以将等式(1)写成
为了求解上述方程式，首先进行离散化处理。将坐标
离散化为 个点，用 来表示，若 ，则有 另外，还需要对方程式进行差分格式处理，二阶导数处理如下
因此，
的二阶导数矩阵可以写成 矩阵矩阵大小为 而不是 ，因为函数的二阶导数无法在终点进行计算，即 和 。动能矩阵为 ，势能矩阵为对角矩阵，即 。则我们现在可以将哈密顿矩阵定义为 。用于生成哈密顿矩阵的代码是
% Make Second Derivative Matrix ---------------
因此，矩阵形式的薛定谔方程是
。MATLAB有内置函数可以求解本征值问题，其代码为
[e_funct, e_values] = eig(H)
其中e_funct是具有对应于第n个本征函数的第n列的
矩阵，并且e_values是按递增顺序的N个本征值的列向量。一般可以求解出N-2个本征值和本征函数，但只有e_values的负值才有意义。为了获得完整的特征向量，我们需要包括端点，其中 。
接下来举一个例子，计算有限深方势阱问题。如图所示是求解得到的有限深方势阱的本征能量谱。
同时还可以求出本征函数以及几率分布，如图所示，
我们还可以改变势函数去计算各种各样的定态束缚态问题，也可以去计算已知解的问题，以验证此方法的可靠性。
程序：
clear;clc;tic;num = 2001; % Number of data points (odd number)% Constants -----------------hbar = 1.055e-34;e = 1.602e-19;m = 9.109e-31;eps0 = 8.854e-12;Ese = 1.6e-19; % Energy scaling factorLse = 1e-9;% Length scaling factorCse = -hbar^2/(2*m) / (Lse^2*Ese); % Schrodinger Eq constant% Potential well parametersU = zeros(num,1);U_matrix = zeros(num-2);% Potential Wells square well% Enter energies in eV and distances in nmxMin = -0.1;xMax = +0.1;x1 = 0.05; % 1/2 well widthU1 = -400; % Depth of well (eV)x = linspace(xMin,xMax, num);for cn = 1 : numif abs(x(cn)) <= x1U(cn) = U1;endends = sprintf('Potential Well: SQUARE');% Graphics -----------------------figure(1);set(gcf,'Name','Potential Energy','NumberTitle','off')plot(x,U,'LineWidth',3);axis([xMin-eps xMax min(U)-50 max(U)+50]);title(s);xlabel('x (nm)');ylabel('energy (eV)');grid on% Make potential energy matrixdx = (x(2)-x(1));dx2 = dx^2;for cn =1:(num-2)U_matrix(cn,cn) = U(cn+1);end% Make Second Derivative Matrixoff = ones(num-3,1);SD_matrix = (-2*eye(num-2) + diag(off,1) + diag(off,-1))/dx2;% Make KE MatrixK_matrix = Cse * SD_matrix; % Make Hamiltonian MatrixH_matrix = K_matrix + U_matrix;% Find Eignevalues E_n and Eigenfunctions psi_N[e_funct, e_values] = eig(H_matrix);% All Eigenvalues 1, 2 , ... n where E_N < 0flag = 0;n = 1;while flag == 0E(n) = e_values(n,n);if E(n) > 0flag = 1;end % ifn = n + 1;end % whileE(n-1) = [];n = n-2;% Corresponding Eigenfunctions 1, 2, ... ,n: Normalizing the wavefunctionfor cn = 1 : npsi(:,cn) = [0; e_funct(:,cn); 0]; area = simpson1d((psi(:,cn) .* psi(:,cn))',xMin,xMax);psi(:,cn) = psi(:,cn)/sqrt(area); % normalizeprob(:,cn) = psi(:,cn) .* psi(:,cn);if psi(5,cn) < 0psi(:,cn) = -psi(:,cn); end % curve starts positiveend % for% Display eigenvalues in Command Windowdisp(' ');disp('========================= ');disp(' ');fprintf('No. bound states found = %0.0g n',n);disp(' ');disp('Quantum State / Eigenvalues En (eV)');for cn = 1 : nfprintf(' %0.0f ',cn);fprintf(' %0.5g n',E(cn));enddisp(' ')disp(' ');% Plot energy spectrumxs(1) = xMin;xs(2) = xMax;figure(2);set(gcf,'Units','Normalized');set(gcf,'Position',[0.5 0.1 0.4 0.6]);set(gcf,'Name','Energy Spectrum','NumberTitle','off')set(gcf,'color',[1 1 1]);set(gca,'fontSize',12);plot(x,U,'b','LineWidth',2);xlabel('position x (nm)','FontSize',12);ylabel('energy U, E_n (eV)','FontSize',12);h_title = title(s);set(h_title,'FontSize',12);hold oncnmax = length(E);for cn = 1 : cnmaxys(1) = E(cn);ys(2) = ys(1);plot(xs,ys,'r','LineWidth',2);end %for axis([xMin-eps xMax min(U)-50 max(U)+50]);% Plots first 5 wavefunctions & probability density functionsif n < 6nMax = n;elsenMax = 5;endfigure(3)clfset(gcf,'Units','Normalized');set(gcf,'Position',[0.05 0.1 0.4 0.6]);set(gcf,'NumberTitle','off');set(gcf,'Name','Eigenvectors & Prob. densities');set(gcf,'Color',[1 1 1]);%nMax = 8;for cn = 1:nMaxsubplot(nMax,2,2*cn-1);y1 = psi(:,cn) ./ (max(psi(:,cn)-min(psi(:,cn))));y2 = 1 + 2 * U ./ (max(U) - min(U));plot(x,y1,'lineWidth',2)hold onplot(x,y2,'r','lineWidth',1)%plotyy(x,psi(:,cn),x,U);axis off%title('psi cn);title_m = ['psi n = ', num2str(cn)] ;title(title_m,'Fontsize',10);subplot(nMax,2,2*cn);y1 = prob(:,cn) ./ max(prob(:,cn));y2 = 1 + 2 * U ./ (max(U) - min(U));plot(x,y1,'lineWidth',2)hold onplot(x,y2,'r','lineWidth',1)title_m = ['psi^2 n = ', num2str(cn)] ;title(title_m,'Fontsize',10);axis offendtoc
m文件simpson1d.m
function integral = simpson1d(f,a,b)% [1D] integration - Simpson's 1/3 rule% f function a = lower bound b = upper bound% Must have odd number of data points% Simpson's coefficients 1 4 2 4 ... 2 4 1numS = length(f); % number of data pointsif mod(numS,2) == 1sc = 2*ones(numS,1);sc(2:2:numS-1) = 4;sc(1) = 1; sc(numS) = 1;h = (b-a)/(numS-1);integral = (h/3) * f * sc;else integral = 'Length of function must be an ODD number' end
参考资料：
http://www.physics.usyd.edu.au/teach_res/mp/mphome.htm
