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

d = 5; % 板の１辺の長さ
f = 20; % 発熱量
h = d/(m - 1); % 分割の幅
n = m^2; % 節点数
c = dt/h^2;
u = 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
      u(i) = 0;
    else % 板の内部
      u(i) = (1 - 4*c)*u0(i) ...
      + c*(u0(i+1) + u0(i-1) + u0(i+m) + u0(i-m)) + dt*f;
    end
  end
end