clear;

// *** 入力パラメータ ***
n = 256;		// 位置の分割数
vp = 10000.0;		// ポテンシャルの高さ
xm = 1.0;
wp = 1 / 32 + %eps;

// *** ハミルトニアンの定義 ***
// エネルギー演算子
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.';
// ポテンシャル
V = zeros(H);
for k = 1:n-1 do
  if ((abs(k * h - 1 / 16) <= wp) | (abs(k * h - 3 / 16) <= wp) |..
      (abs(k * h - 5 / 16) <= wp) | (abs(k * h - 7 / 16) <= wp) | ..
      (abs(k * h - 9 / 16) <= wp) | (abs(k * h - 11 / 16) <= wp) | ..
      (abs(k * h - 13 / 16) <= wp) | (abs(k * h - 15 / 16) <= wp)) then
    V(k,k) = vp;
  end
end
// 全ハミルトニアン
H = H + V;

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

// *** グラフのプロット ***
x = linspace(h, xm - h, n - 1);
// ポテンシャルのプロット
xsetech([0,0,1,0.5]);
plot(x,diag(V),'--m');
zoom_rect([0,0,xm,vp*1.2]);
xlabel("Position");
ylabel("Potential");
// 波動関数のプロット
xsetech([0,0.475,1,0.5]);
plot(x,Phi(:,1),'-r');
plot(x,Phi(:,2),'-g');
plot(x,Phi(:,3),'-b');
xlabel("Position");
ylabel("Wave function");