clear;

vw = 0.064;			// ポテンシャルの幅
vh = - (70.7 * %pi) ^ 2;	// ポテンシャルの深さ
dim = 500;			// 空間の分割数
h = 2.0 / dim;			// 空間の刻み幅

// *** 位置ベクトル ***
X = linspace(-1,1,dim + 1);

// *** ハミルトニアンの定義 ***
// エネルギー演算子
vd = 2 / h ^ 2;
vdt = - 1 / h ^ 2;
SD = diag(vdt * ones(1,dim),1);
D = diag(vd * ones(1,dim + 1));
H = SD + D + SD.';
// ポテンシャル
V = zeros(H);
for i = 1:dim + 1 do
  if abs(X(i)) < vw / 2 then
    V(i,i) = vh;
  end
end
// 全ハミルトニアン
H = sparse(H + V);		// 疎行列へ変換するほうが実行速度が速い

// *** 常微分方程式の定義 ***
function dphi = schiff(t,phi)
  dphi = - %i * H * phi;
endfunction

// *** 波束の初期設定 ***
xe = - 0.3;			// 初期波束の位置の期待値(波束の中心位置)
dx = 0.035;			// 初期波束の空間的広がり
pe = 50 * %pi;			// 初期波束の運動量の期待値
phi0 = zeros(dim + 1,1);
for k = 1:dim + 1 do
  // 初期波束
  phi0(k) = exp(-(X(k) - xe) ^ 2 / (4 * dx ^ 2) + %i * pe * X(k)) ..
           / (2 * %pi * (dx) ^ 2) ^ 0.25;
end

// *** 時間刻みの設定 ***
dt = h ^ 2 / 4;
nstep = 700;
itvl = 10;
tf = dt * itvl;
T = [0:dt:tf];

// *** プロットと画像出力 ***
driver('GIF');				// GIF形式で出力

// 初期値のプロット
xinit('schiff000.GIF');			// 画像ファイルの指定
plot(X,diag(V));			// ポテンシャルのプロット
plot(X, real(phi0),'-g');		// 波動関数の実部のプロット
plot(X, abs(phi0),'-r');		// 波動関数の絶対値のプロット
// ステップ数を表示
xstring(0.3,4,"step = 0");
// 目盛りを非表示にする
g = gca();
g.axes_visible = 'off';
zoom_rect([-0.5,-5,0.5,5]);
xend();					// 画像ファイルを閉じる

for k = 1:nstep/itvl do
  // ファイル番号の先頭にゼロを補う
  filenum = string(k);
    if length(filenum) == 1 then
    filenum = strcat(["00",filenum]);
  elseif length(filenum) == 2 then
    filenum = strcat(["0",filenum]);
  end
  // 画像ファイルの指定
  xinit(strcat(['schiff',filenum,'.GIF']));
  // 微分方程式の数値解
  phi = ode(phi0,0,T,schiff);
  // ポテンシャルと波動関数のプロット
  plot(X,diag(V));			// ポテンシャルのプロット
  plot(X,real(phi(:,itvl)),'-g');	// 波動関数の実部のプロット
  plot(X,abs(phi(:,itvl)),'-r');	// 波動関数の絶対値のプロット
  // ステップ数を表示
  stepnum = string(k * itvl);
  xstring(0.3,4,strcat(["step = ",stepnum]));
  // 目盛りを非表示にする
  g = gca();
  g.axes_visible = 'off';
  zoom_rect([-0.5,-5,0.5,5]);
  xend();				// 画像ファイルを閉じる
  phi0 = phi(:,itvl);
end

// グラフ出力先をコンピュータ画面出力へもどす
driver('X11');