% 長方形領域における熱方程式 u_t=△u (Dirichlet境界条件) を解くための差分方程式
%   A U^{n+1}=B U^n
% の行列 A, B を求める。
% Nx=3; Ny=3; hx=1/Nx; hy=1/Ny; theta=0.5; tau=0.5/(1/hx^2+1/hy^2); lamx=tau/hx^2; lamy=tau/hy^2;
% heat2n_mat(Nx,Ny,lamx,lamy,theta)
% A =
%    1.5000   -0.1250   -0.1250         0
%   -0.1250    1.5000         0   -0.1250
%   -0.1250         0    1.5000   -0.1250
%         0   -0.1250   -0.1250    1.5000
% B =
%    0.5000    0.1250    0.1250         0
%    0.1250    0.5000         0    0.1250
%    0.1250         0    0.5000    0.1250
%         0    0.1250    0.1250    0.5000
function [A,B]=heat2d_mat(Nx,Ny,lambdax,lambday,theta)
  Ix=speye(Nx+1,Nx+1);
  Iy=speye(Ny+1,Ny+1);
  vx=[sqrt(2); ones(Nx-2,1); sqrt(2)];
  Jx=sparse(diag(vx,1)+diag(vx,-1));
  vy=[sqrt(2); ones(Ny-2,1); sqrt(2)];
  Jy=sparse(diag(vy,1)+diag(vy,-1));
  Kx=2*Ix-Jx;
  Ky=2*Iy-Jy;
% column first
  A=kron(Ix,Iy)+theta*lambday*kron(Ix,Ky)+theta*lambdax*kron(Kx,Iy);
  B=kron(Ix,Iy)-(1-theta)*lambday*kron(Ix,Ky)-(1-theta)*lambdax*kron(Kx,Iy);
% row first
%  A=kron(Iy,Ix)+theta*lambdax*kron(Iy,Kx)+theta*lambday*kron(Ky,Ix);
%  B=kron(Iy,Ix)-(1-theta)*lambdax*kron(Iy,Kx)-(1-theta)*lambday*kron(Ky,Ix);
  
