clear;

// *** 入力パラメータ ***
n = 256;		// 位置の分割数
v0 = 100.0;		// ポテンシャルの高さ
xm = 1.0;		// 位置の最大値
w0 = xm / 16;		// 中央ポテンシャル壁の幅

// *** ハミルトニアンの定義 ***
// エネルギー演算子
h = xm / n;		// 位置の刻み幅
vd = 2 / h / h;
vdt = - 1 / h / h;
SD = diag(vdt * ones(1,n - 2),1);
D = diag(vd * ones(1,n - 1));
H = SD + D + SD.';
// ポテンシャル
for i = 1:n - 1 do
  if abs(i * h - xm / 2) <= w0 then
  H(i,i) = H(i,i) + v0;
  end
end

// *** 固有値問題を解く ***
// Phi: 固有ベクトル(波動関数)
// E: 固有値(固有エネルギー)の対角行列
[Phi,E] = spec(H);

// *** グラフのプロット ***
x = linspace(xm / n, xm - xm / n, n - 1);
// 波動関数のプロット
plot(x,Phi(:,1),'-+r');
plot(x,Phi(:,2),'-xg');
plot(x,Phi(:,3),'-*b');
// ポテンシャルのプロット
ptx1 = xm / 2 - w0;			// 中央壁の端
ptx2 = xm / 2 + w0;			// 中央壁の端
ptv = v0 / 2000;			// ポテンシャルの高さ適当に規格化
ptx = [0, ptx1, ptx1, ptx2, ptx2, xm];
pt = [0, 0, ptv, ptv, 0, 0];
plot(ptx,pt,'--m');			// ポテンシャル壁のプロット
legend(['E1','E2','E3','Potential'],2);