function u = theat2d3(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,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.0;
      b(i) = 0;
    else % 板の内部
       A(i,i) = 1 + 2*c; % u(i)の係数
       A(i,i-1) = -0.5*c; % u(i-1)の係数
       A(i,i+1) = -0.5*c; % u(i+1)の係数
       A(i,i-m) = -0.5*c; % u(i-m)の係数
       A(i,i+m) = -0.5*c; % u(i+m)の係数
       b(i) = (1 - 2*c)*u0(i) ...
       + 0.5*c*(u0(i+1) + u0(i-1) + u0(i+m) + u0(i-m)) + dt*f;
    end
  end
end
u = A\b;