我正在尝试数值求解Klein-Gordon方程，它可以找到here。为了确保我正确地解决了它，我将它与可以在同一链接上找到的解析解进行比较。我使用的是有限差分法和MatLab。初始空间条件是已知的，而不是初始时间条件。
我首先初始化常量和时空坐标系：

close all
clear 
clc
 
%% Constant parameters

A      = 2; 
B      = 3; 
lambda = 2; 
mu     = 3; 
a      = 4; 
b      = - (lambda^2 / a^2) + mu^2;  

%% Coordinate system

number_of_discrete_time_steps = 300; 
t = linspace(0, 2, number_of_discrete_time_steps); 
dt = t(2) - t(1); 

number_of_discrete_space_steps = 100; 
x = transpose( linspace(0, 1, number_of_discrete_space_steps) ); 
dx = x(2) - x(1);

接下来，我定义并绘制分析解决方案：

%% Analitical solution

Wa = cos(lambda * x) * ( A * cos(mu * t) + B * sin(mu * t) );

figure('Name', 'Analitical solution'); 
surface(t, x, Wa, 'edgecolor', 'none'); 
colormap(jet(256)); 
colorbar; 
xlabel('t'); 
ylabel('x'); 
title('Wa(x, t) - analitical solution');

解析解的曲线图如here所示。
最后，我定义了初始空间条件，执行有限差分法算法并绘制了解：

%% Numerical solution

Wn = zeros(number_of_discrete_space_steps, number_of_discrete_time_steps);

Wn(1, :) = Wa(1, :); 
Wn(2, :) = Wa(2, :); 

for j = 2 : (number_of_discrete_time_steps - 1)  

    for i = 2 : (number_of_discrete_space_steps - 1) 

        Wn(i + 1, j) = dx^2 / a^2 ...                        
                     * ( ( Wn(i, j + 1) - 2 * Wn(i, j) + Wn(i, j - 1) ) / dt^2 + b * Wn(i - 1, j - 1) ) ...
                     + 2 * Wn(i, j) - Wn(i - 1, j); 
    end

end

figure('Name', 'Numerical solution'); 
surface(t, x, Wn, 'edgecolor', 'none'); 
colormap(jet(256)); 
colorbar; 
xlabel('t'); 
ylabel('x'); 
title('Wn(x, t) - numerical solution');

数值解的曲线图如here所示。
绘制的两个曲线图不一样，这证明我在算法中做错了什么。问题是，我找不到错误。请帮我找到他们。
总而言之，请帮助我更改代码，以便两个绘制的图形变得大致相同。谢谢您抽时间见我。