Solving coupled pdes with FFT

I have been trying to generalise some working code from Professor Steve Brunton which solves the one way wave pde using the FFT method. This is from his book at http://databookuw.com The explanation of the method and the code is at https://www.youtube.com/watch?v=RilggJd1_LY I paste the code below and it works very well.

I would like to solve a pair of coupled 'two-way' wave equations by the same method. In fourier space the equivalent first order system is

uhat_t = phat

vhat_t = qhat

phat_t = mu^2 k^2 uhat + qhat

qhat_t - lambda^2 mu^2 k^2 vhat -phat

where lambda and mu are real constants, k is the spatial frequency and the fourier transformed vector of unknowns is [uhat,vhat,phat,qhat]

Clearly, there is a vector of unknowns instead of just one scalar unknown in the working example. It is not clear to me that ode45 can be used for this problem as the 2D array of output for the scalar problem would somehow need to be extended to a 3D array with the use of a 'page' dimension. I have attemepted this but cannot get anywhere.

Does anyone have an example of such a higher dimensional problem or is it possible to adapt the code below?

Many thanks in anticipation of your help.

clear all, close all, clc

% Define spatial domain

c = 2; % Wave speed

L = 20; % Length of domain

N = 1000; % Number of discretization points

dx = L/N;

x = -L/2:dx:L/2-dx; % Define x domain

% Define discrete wavenumbers

kappa = (2*pi/L)*[-N/2:N/2-1];

kappa = fftshift(kappa'); % Re-order fft wavenumbers

% Initial condition

u0 = sech(x);

uhat0 = fft(u0);

% Simulate in Fourier frequency domain

dt = 0.025;

t = 0:dt:100*dt;

[t,uhat] = ode45(@(t,uhat)rhsWave(t,uhat,kappa,c),t,uhat0);

% Alternatively, simulate in spatial domain

[t,u] = ode45(@(t,u)rhsWaveSpatial(t,u,kappa,c),t,u0);

% Inverse FFT to bring back to spatial domain

for k = 1:length(t)

u(k,:) = ifft(uhat(k,:));

end

% Plot solution in time

subplot(1,2,1)

h=waterfall(real(u(1:10:end,:)));

set(h,'LineWidth',2,'FaceAlpha',.5);

colormap(jet/1.5)

subplot(1,2,2)

imagesc(flipud(real(u)));

colormap jet

%% FIGURES (PRODUCTION)

figure

CC = colormap(jet(101));

dt = 0.025;

for k = 1:length(t)

u(k,:) = real(ifft(uhat(k,:)));

if(mod(k-1,10)==0)

plot(x,u(k,:),'Color',CC(k,:),'LineWidth',1.5)

hold on, grid on, drawnow

end

end

% xlabel('Spatial variable, x')

% ylabel('Temperature, u(x,t)')

axis([-L/2 L/2 -.1 1.1])

set(gca,'LineWidth',1.2,'FontSize',12);

set(gcf,'Position',[100 100 550 220]);

set(gcf,'PaperPositionMode','auto')

%

figure

subplot(1,2,1)

h=waterfall(u(1:10:end,:));

set(h,'LineWidth',2,'FaceAlpha',.5);

colormap(jet/1.5)

% view(22,29)

view(3,21)

set(gca,'LineWidth',1.5)

set(gca,'XTick',[0 500 1000],'XTickLabels',{})

set(gca,'ZTick',[0 .5 1],'ZTickLabels',{})

set(gca,'YTick',[0 5 10],'YTickLabels',{})

subplot(1,2,2)

imagesc(flipud(u));

set(gca,'LineWidth',1.5)

set(gca,'XTick',[0 500 1000],'XTickLabels',{})

set(gca,'YTick',[0 50 100],'YTickLabels',{})

colormap jet

set(gcf,'Position',[100 100 600 250])

set(gcf,'PaperPositionMode','auto')

__________________________________________________

function duhatdt = rhsWave(t,uhat,kappa,c)

duhatdt = -c*i*kappa.*uhat;

________________________________________________

function dudt = rhsWaveSpatial(t,u,kappa,c)

uhat = fft(u);

duhat = i*kappa.*uhat;

du = ifft(duhat);

dudt = -c*du;