function u = theat2d2(u0,m,dt)
% 熱伝導シミュレーション（２次元、非定常、陰解法）
%  m: x,y方向の分割数, u0: 初期値, dt: 時間の刻み幅
%  u: １ステップ後の温度

d = 5; % 板の１辺の長さ
f = 20; % 発熱量
h = d/(m - 1); % 分割の幅
n = m^2; % 節点数
c = dt/h^2;
A = zeros(n); b = zeros(n,1); % 初期化

for jj=1:m
  for ii=1:m
    i = ii + (jj - 1)*m;
    if (ii == 1) || (ii == m) || (jj == 1) || (jj == m)
       % 境界条件 u = 0
       A(i,i) = 1; b(i) = 0;
    else % 板の内部
       A(i,i) = 1 + 4*c; % u(i)の係数
       A(i,i-1) = -c; % u(i-1)の係数
       A(i,i+1) = -c; % u(i+1)の係数
       A(i,i-m) = -c; % u(i-m)の係数
       A(i,i+m) = -c; % u(i+m)の係数
       b(i) = u0(i) + dt*f; % 右辺
    end
  end
end
u = A\b; % Au=bを解く