我用我认为是 Crank-Nicolson 算法的方式实现了 2D 热方程，方法如下：

function[newDomain] = heatCN(domain, dt, Nx, Ny, hx, hy, spanX, spanY)
aux = domain;
newDomain = domain;
for i=2:Nx-1
    for j=2:Ny-1
        derX1 = dt * (-2.0 * domain(i,j) + domain(i+1, j) + domain(i-1,j)) / hx^2;
        derY1 = dt * (-2.0 * domain(i,j) + domain(i, j+1) + domain(i, j-1)) / hy^2;
        aux(i,j) = derX1 + derY1 + domain(i,j);
    end;
end;

for i=2:Nx-1
    for j=2:Ny-1
        derX2 = dt * (-2.0 * domain(i,j) + domain(i+1, j) + domain(i-1,j)) / hx^2;
        derY2 = dt * (-2.0 * domain(i,j) + domain(i, j+1) + domain(i, j-1)) / hy^2;
        derXT = dt * (-2.0 * aux(i,j) + aux(i+1, j) + aux(i-1,j)) / hx^2;
        derYT = dt * (-2.0 * aux(i,j) + aux(i, j+1) + aux(i, j-1)) / hy^2;
        newDomain(i,j) = domain(i,j) + 0.5 * (derX2 + derY2 + derXT + derYT);
    end;
end;
end;

据我所知，CN 方法应该是无条件稳定的，这意味着对于任何（体面的）时间步，它都应该给我正确的输出。但是我在这个程序中看到的行为是不稳定的。