举一个matlab绘制三维反应扩散方程图灵斑图的详细例子，要求用差分法表示反应扩散方程

在MATLAB中，绘制三维反应扩散方程的图灵斑图案通常涉及离散化连续偏微分方程，例如欧拉方法或有限差分法。以下是一个使用有限差分法简化版本的示例，展示如何创建一个简单的图灵斑模式：

% 定义函数参数
L = 10; % 空间域范围
dx = dy = 0.1; % 步长
dt = 0.1; % 时间步长
Nt = 500; % 时间迭代次数
D = 1; % 扩散系数
kappa = 0.2; % 图灵稳定性参数
r = 0.5; % 初始波幅
theta = 0; % 初始相位

% 初始化网格和状态矩阵
[x, y] = meshgrid(0:L:dx, 0:L:dy);
u = zeros(size(x));
u(x.^2 + y.^2 <= r^2) = 1; % 圆形区域作为初始波

% 用二维数组表示空间和时间
[u_t, x_t, y_t] = ndgrid(linspace(0, L, length(x)+1)', ... 
                            linspace(0, L, length(y)+1)', 
                            linspace(0, Nt, dt*Nt+1));

% 有限差分计算反应扩散方程
for t = 1:Nt
    % 计算扩散项
    du_dx = (u_t(:,:,t+1) - u_t(:,:,t)) ./ dx;
    du_dy = (u_t(:,end,:)-u_t(:,1:end-1,:))./dy;
    
    % 计算图灵稳定性项
    term1 = kappa * (1 - 2*u_t(:,:,t).^2);
    term2 = sin(theta + 2*pi*(x_t + y_t));
    
    % 更新状态矩阵
    u_t(:,:,t+1) = u_t(:,:,t) + D * (du_dx + du_dy) + term1 .* term2;
end

% 将结果转换为表面图
figure
surf(x_t(2:end-1,2:end-1,:), y_t(2:end-1,2:end-1,:), u_t(2:end-1,2:end-1,:))
view(2); % 显示侧面视角
title('反应扩散方程的图灵斑')
xlabel('x')
ylabel('y')
zlabel('浓度')