clear;

// *** 入力パラメータ ***
n = 256;		// 位置の分割数
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));
A = SD + D + SD.';
// ポテンシャル
function H = potential(vp)
  // ポテンシャル
  V = zeros(A);
  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 = A + V;
endfunction

// *** 波動関数・固有エネルギーの計算とプロット ***
// VP = 50
H = potential(50);
[Phi,E] = spec(H);
plot(diag(E),'-+r');
// VP = 250
H = potential(250);
[Phi,E] = spec(H);
plot(diag(E),'-xg');
// VP = 1250
H = potential(1250);
[Phi,E] = spec(H);
plot(diag(E),'-*b');
// VP = 6250
H = potential(6250);
[Phi,E] = spec(H);
plot(diag(E),'-sm');
// VP = 31250
H = potential(31250);
[Phi,E] = spec(H);
plot(diag(E),'-sc');

// *** グラフの装飾 ***
a = gca();
a.data_bounds = [0,1E1; 50,1E5];	// 表示範囲の設定
a.log_flags = "nl";			// 縦軸をログスケールに
// 凡例
legend(['VP=50.0','VP=250.0','VP=1250.0','VP=6250.0','VP=31250.0'],4);
// 軸ラベル
xlabel("Eigen States");
ylabel("Ek");