clear;

// *** “º‚Ì•¨�«’l ***
k = 398;		// ”M“`“±“x 398 W/m/K
rho = 8.92E3;		// –§“x 8.92�~10^3 kg/m^3

cpm = 24.44;		// ’è�󂃋”ä”M 24.44 J/mol/K
m = 63.546;		// Œ´Žq—Ê g/mol
mkg = m * 1E-3;		// Œ´Žq—Ê kg/mol
cp = cpm / mkg;		// ’è�ÏŽ¿—Ê”ä”M J/kg/K

d2 = k / (rho * cp)	// ”MŠgŽU—¦ m^2/s

// *** –³ŽŸŒ³‰»‘O‚ÌŽžŠÔ‚Æ‹óŠÔ ***
lambda = 1.5;		// –_‚Ì’·‚³ ƒÉ (m)
tdmax = 60 * 60;	// t'max (s) ‚Ü‚ÅŒvŽZ

// *** –³ŽŸŒ³‰» ***
tau = lambda ^ 2 / d2;	// ŽžŠÔ‚Ì–³ŽŸŒ³‰»ŒW�”
tmax = tdmax / tau;	// –³ŽŸŒ³‰»ŽžŠÔ‚Å tmax ‚Ü‚ÅŒvŽZ
temp0 = 100;		// u=1 ‚ɑ΂µ‚ÄT'=100�Ž

// *** ŽžŠÔ‚Æ‹óŠÔ‚Ì•ªŠ„ ***
n = 50;		// ‹óŠÔ‚Ì•ªŠ„�”
dx = 1 / (n - 1);	// ‹óŠÔ‚Ì��‚Ý•�

r = 1 / 6;

dt = r * dx ^ 2;	// ŽžŠÔ‚Ì��‚Ý•�
m = round(tmax / dt);	// ŽžŠÔ‚Ì•ªŠ„�”

// *** �‰Šú�ðŒ�‚Æ‹«ŠE�ðŒ� ***
// �‰Šú�ðŒ�
U = zeros(n,m);
// ‹«ŠE�ðŒ�
U(1,1) = 0;
U(n,1) = 1;

// *** �s—ñP ***
s = 1 - 2 * r;
P = toeplitz([s, r, zeros(1, n - 2)]);

// *** •Î”÷•ª•û’öŽ®‚ÌŒvŽZ ***
for j = 1:(m - 1)
  // ‹«ŠEˆÈŠO‚ÌŒvŽZ
  U(:, j + 1) = P * U(:,j);
  // ‹«ŠE�ðŒ�
  U(1, j + 1) = 0;
  U(n, j + 1) = 1;
end

// *** ƒOƒ‰ƒt‚̃vƒ�ƒbƒg ***
X = linspace(0,lambda,n);
plot(X,temp0 * U(:,m),'o-b');
zoom_rect([0,0,lambda,temp0]);

xlabel("Position (m)");
ylabel("Temperature (degree C)");