给我1个高精度隐式有限差分法求解二维非线性Klein-Gordon方程的matlab代码

以下是一个基于高精度隐式有限差分法求解二维非线性Klein-Gordon方程的MATLAB代码。请注意，此代码仅供参考，可能需要根据具体问题进行修改。

clear;clc;
%定义参数
L = 50; %空间区域长度
h = 0.1; %网格尺寸
N = L/h; %网格数
dt = 0.01; %时间步长
t_max = 100; %模拟时间
x = linspace(-L/2,L/2,N); %空间网格
y = linspace(-L/2,L/2,N); %空间网格
[X,Y] = meshgrid(x,y); %生成网格
k = 1; %常数
mu = 1; %常数

%定义初始值
u = exp(-(X.^2+Y.^2)/2); %初始条件
v = zeros(size(u)); %初始速度

%构造矩阵
e = ones(N,1);
D = spdiags([-e 2*e -e], [-1 0 1], N, N);
D(1,end) = -1;
D(end,1) = -1;
D = D/h^2;

%时间演化
for n = 1:round(t_max/dt)
    u_old = u;
    v_old = v;
    
    %求解非线性项
    u2 = u_old.^2;
    v2 = v_old.^2;
    uv = u_old.*v_old;
    f = -k^2*(1+mu*u2).*u_old + v_old;
    g = -k^2*(1+mu*v2).*v_old + u_old;
    F = [f(:);g(:)];
    
    %构造矩阵
    A = [speye(N^2)+dt*D, sparse(N^2,N^2);sparse(N^2,N^2), speye(N^2)+dt*D];
    B = [speye(N^2)-dt*D, sparse(N^2,N^2);sparse(N^2,N^2), speye(N^2)-dt*D];
    
    %求解线性方程组
    b = B*[u_old(:);v_old(:)] + dt*F;
    x = A\b;
    
    %更新解
    u = reshape(x(1:N^2),[N,N]);
    v = reshape(x(N^2+1:end),[N,N]);
    
    %绘图
    surf(X,Y,u,'edgecolor','none');
    axis([-L/2 L/2 -L/2 L/2 -1 1]);
    title(['t=',num2str(n*dt)]);
    drawnow;
end