编写一段用有限差分法求解对流扩散方程的MATLAB代码，并将结果可视化

在MATLAB中，我们可以使用有限差分方法（Finite Difference Method, FDM）来解决二维对流扩散方程。这里是一个简单的例子，我们将考虑二维空间和时间的一维平滑对流传导方程。假设我们要在单元网格上计算：

∂u/∂t = D * (∂²u/∂x² + ∂²u/∂y²) - u

其中，u(x,y,t)是浓度分布，D是扩散系数。

% 定义常量和变量
L = 1; % 空间域长度
N = 50; % 网格点数
T = 1; % 时间总长
dt = T / N; % 时间步长
dx = L / (N - 1); % 空间步长
D = 1; % 扩散系数

% 初始化网格和初始条件
[x, y] = meshgrid(0:dx:L, 0:dy:L);
u0 = sin(pi*x).*sin(pi*y); % 可以替换为任意初始分布

% 创建储存结果的矩阵
[u, t] = deal(zeros(N, N, N));

% 使用向前欧拉法和中心差异法计算
for n = 1:N
    u(:,:,n) = u0;
    for i = 2:N-1
        for j = 2:N-1
            u(i,j,n+1) = u(i,j,n) + dt/D * (-u(i+1,j,n) - u(i-1,j,n) - u(i,j+1,n) - u(i,j-1,n) + 4*u(i,j,n));
        end
    end
end

% 结果可视化
surf(u(:, :, int2str(end))) % 绘制最终时刻的浓度分布
xlabel('X')
ylabel('Y')
zlabel('Concentration')
title('Solution at Time T')

这段代码首先初始化了一个二维网格和初始浓度分布，然后使用迭代的方式逐步更新每个网格点的浓度值。最后，`surf`函数用于三维可视化结果。