function u = theat1d3(u0,n,dt)
% 熱伝導シミュレーション（１次元、非定常、クランク・ニコルソン法）
%  u0: 温度の初期値, n: 節点の数, dt: 時間の刻み幅

d = 5; % 棒の長さ
f = 20; % 発熱量
h = d/(n - 1); % 分割の幅
c = dt/h^2;
A = zeros(n); b = zeros(n,1); % 初期化
A(1,1) = 1; b(1) = 0; % u(1) = 0
A(n,n) = 1; b(n) = 0; % u(n) = 0
for i=2:n-1
  A(i,i) = 2 + 2*c;
  A(i,i-1) = -c;
  A(i,i+1) = -c;
  b(i) = (2 - 2*c)*u0(i) + c*u0(i-1) + c*u0(i+1) + 2*dt*f;
end
u = A\b;