给我一个有限差分法求解二维非线性Klein-Gordon方程的matlab求解代码

以下是一个基于有限差分法求解二维非线性Klein-Gordon方程的Matlab求解代码：

% 设置求解区域和时间范围
Lx = 50; % 求解区域的长度
Ly = 50; % 求解区域的宽度
T = 100; % 求解时间的范围

% 设置网格参数
Nx = 101; % x方向上的网格数
Ny = 101; % y方向上的网格数
Nt = 1000; % 时间步数
dx = Lx/(Nx-1); % x方向上的网格间距
dy = Ly/(Ny-1); % y方向上的网格间距
dt = T/Nt; % 时间步长

% 初始化变量
u = zeros(Nx,Ny,Nt); % 存储解
x = linspace(0,Lx,Nx);
y = linspace(0,Ly,Ny);
t = linspace(0,T,Nt);

% 设置初始条件
[X,Y] = meshgrid(x,y);
u(:,:,1) = exp(-((X-Lx/2).^2+(Y-Ly/2).^2)/10);

% 设置边界条件
u(1,:,:) = u(Nx,:,:);
u(Nx,:,:) = u(1,:,:);
u(:,1,:) = u(:,Ny,:);
u(:,Ny,:) = u(:,1,:);

% 求解
for n = 1:Nt-1
    % 使用中心差分求解x方向上的二阶导数
    uxx = (u(1:Nx-2,:,n)-2*u(2:Nx-1,:,n)+u(3:Nx,:,n))/(dx^2);
    % 使用中心差分求解y方向上的二阶导数
    uyy = (u(:,1:Ny-2,n)-2*u(:,2:Ny-1,n)+u(:,3:Ny,n))/(dy^2);
    % 计算非线性项
    f = u(:,:,n).^3-u(:,:,n);
    % 更新解
    u(2:Nx-1,2:Ny-1,n+1) = u(2:Nx-1,2:Ny-1,n)+dt*(uxx+uyy+f(2:Nx-1,2:Ny-1));
    % 使用边界条件更新解
    u(1,:,:) = u(Nx,:,:);
    u(Nx,:,:) = u(1,:,:);
    u(:,1,:) = u(:,Ny,:);
    u(:,Ny,:) = u(:,1,:);
end

% 绘制解在不同时间的图像
for n = 1:100:Nt
    contourf(x,y,u(:,:,n),20);
    colorbar;
    title(['t = ' num2str(t(n))]);
    axis([0 Lx 0 Ly]);
    xlabel('x');
    ylabel('y');
    pause(0.01);
end

在这个代码中，我们使用了中心差分法来求解二维非线性Klein-Gordon方程。在每个时间步长内，我们首先使用中心差分法求解x方向和y方向上的二阶导数，然后计算非线性项，并使用欧拉方法更新解。最后，我们使用边界条件更新解，并绘制解在不同时间的图像。