Parabolic PDE's such as those in the book can usually be solved using the Method of Lines. First you create some mesh for the $x$ direciton. I will assume that you used some uniform spacing since the plots don't show any characteristics that show the need of non-uniformity. Next you recast your equations with only the time derviative on the left hand side and change the derivatives in $x$ to finite difference approximations. Here is the first equation for a general interior point:
$
\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}
$
I am going to leave the second equation, boundary equations, and points near the boundary to you. Now you have a coupled set of ODE's for $i = 1 ... N$. From your initial conditions you can assign $u_{1,i}$ and $u_{2,i}$ at the first time step. Now at each time step you satisfy the above discretized equations at different times based on which time integration algorithm you are using. If it is explicit Euler, you satisfy it at the beginning of each time step. If it is implicit Euler, the end.
In matlab, however, there is an easy way to handle all of these (and many more complicated) methods. What you want is a function that returns a vector of values equal to the right hand side of the above equation given a vector of $u_{j,i}$. If you assume periodic boundary conditions you get:
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
Now you can feed this to any of matlabs built-in ODE solvers. I found that ode15s performed ratherly well. I also assumed sinusoidal ICS, but it doesn't appear to matter.
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);
The results give: 