clear;

// *** 時間と空間の分割 ***
n = 33;			// 空間の分割数
dx = 1 / (n - 1);	// 空間の刻み幅

r = 1 / 6;

dt = r * dx ^ 2;	// 時間の刻み幅
m = 1601;		// 時間の分割数
dm = 200;		// プロットする m の間隔

// *** 初期条件と境界条件 ***
// 初期条件
U = zeros(n,m);
// 境界条件
U(1,1) = 0;
U(n,1) = 1;

// *** 行列P ***
s = 1 - 2 * r;
P = toeplitz([s, r, zeros(1, n - 2)]);

// *** 偏微分方程式の計算 ***
for j = 1:(m - 1)
  // 境界以外の計算
  U(:, j + 1) = P * U(:,j);
  // 境界条件
  U(1, j + 1) = 0;
  U(n, j + 1) = 1;
end

// *** グラフのプロット ***
plot([0:dx:1],U(:,1:dm:m));