书中提到的抛物线偏微分方程通常可以使用直线法求解。首先，您为$x$方向。我将假设您使用了一些均匀的间距，因为这些图没有显示任何表明需要不均匀性的特征。接下来，您仅使用左侧的时间导数重新构建方程，并更改导数$x$到有限差分近似。这是一般内点的第一个方程：

$\frac{\partial u_{1,i}}{\partial t} = \varepsilon (u_{2,i} - u_{1,i}) - \frac{u_{1,i+1} - 2 u_{1,i} + u_{1,i-1}}{\Delta x^2} - \frac{u_{1,i+2} - 4 u_{1,i+1} + 6 u_{1,i} - 4 u_{1,i-1} + u_{1,i-2}}{\Delta x^4} - u_{1,i} \frac{u_{1,i+1} - u_{1,i-1}}{2 \Delta x}$

我将把第二个方程、边界方程和边界附近的点留给你。现在你有一组耦合的 ODE$i = 1 ... N$. 根据您的初始条件，您可以分配$u_{1,i}$和$u_{2,i}$在第一个时间步。现在，在每个时间步，根据您使用的时间积分算法，您在不同时间满足上述离散方程。如果它是显式欧拉，则在每个时间步开始时满足它。如果是隐式欧拉，结束。

然而，在 matlab 中，有一种简单的方法可以处理所有这些（以及许多更复杂的）方法。你想要的是一个函数，它返回一个等于上述等式右侧的值向量，给定一个向量$u_{j,i}$. 如果你假设周期性边界条件，你会得到：

function [ u_prime ] = derivative( t, u, delta_x )

u_prime = zeros(length(u),1);
u = [u(end-3:end); u; u(1:4)];

if t < 200
    epsilon = 0;
else
    epsilon = 0.1;
end


for i = 5:2:length(u) - 4;
    u_prime(i-4) = epsilon*(u(i+1) - u(i)) - ...
        (u(i+2) - 2*u(i) + u(i-2))/delta_x^2 - ...
        (u(i+4) - 4*u(i+2) + 6*u(i) - 4*u(i-2) + u(i-4))/delta_x^4 - ...
        u(i)*(u(i+2) - u(i-2))/(2*delta_x);
    j = i+1;
    u_prime(j-4) = epsilon*(u(j-1) - u(j)) - ...
        (u(j+2) - 2*u(j) + u(j-2))/delta_x^2 - ...
        (u(j+4) - 4*u(j+2) + 6*u(j) - 4*u(j-2) + u(j-4))/delta_x^4 - ...
        u(j)*(u(j+2) - u(j-2))/(2*delta_x); 
    end
end

现在您可以将它提供给任何 matlabs 内置的 ODE 求解器。我发现 ode15s 的表现相当不错。我还假设是正弦 ICS，但这似乎并不重要。

N = 1000;   % Number of space discretizations

x = linspace(0, 150, N);
u_0 = zeros(2*N,1);
u_0(1:2:end-1) = sin(2*x/10);   % u_1
u_0(2:2:end) = -sin(4*x/10);    % u_2
delta_x = x(2) - x(1);

[t,u] = ode15s(@(t,u) derivative(t,u,delta_x), [0 400], u_0);

结果给出：