I have this funny equation
$$ \frac{\partial^2 u}{\partial t^2} = \frac{\partial^3 u}{\partial x^3}, \qquad x \in [0,1], \qquad t \in (0,T] $$ with initial conditions $u(x,0) = \sin(2\pi x)$, $\frac{\partial u}{\partial t}(x,0) = 0$ and periodic boundary conditions: $u(0,t) = u(1,t)$, $\frac{\partial u}{\partial x}(0,t) = \frac{\partial u}{\partial x}(1,t)$ and $\frac{\partial^2 u}{\partial x^2}(0,t) = \frac{\partial^2 u}{\partial x^2}(1,t)$. The exact solution of this model is $$ u(x,t) = \frac{1}{2}[e^{2 \pi^{\frac32}t} \sin (2\pi x - 2\pi^{\frac32} t) + e^{-2\pi^{\frac32}t} \sin (2\pi x + 2 \pi^{\frac32}t)] $$ I would like to use boundary value technique to numerically solve this equation, i.e. i would like to write the following system of $2MN$ unknowns in matrix form and solve it in Matlab.
I partition time as $0 = t_0 < t_1 < \dots < t_{N} = 0.1$. Step size is of length $h=\Delta t =\frac{0.1}{N}$, i.e. $t_n = t_0 + nh$, $n=1,2,\dots,N$.
I partition space as $0 = x_0 < x_1 < \dots < x_{M} = 1$. Step size is of length $h'=\Delta x =\frac{1}{M}$, i.e. $x_i = x_0 + ih'$, $i=1,2,\dots,M$.
I rewrite the equation as system,
\begin{align*} u_t &= v, \\ v_t &= u_{xxx} \end{align*}
and using central approximation for time derivative and second order scheme for space derivative
\begin{align*} \frac{u_i^{n+1} - u_i^{n-1}}{2\Delta t} &= v_i^n \\ \frac{v_i^{n+1} - v_i^{n-1}}{2\Delta t} &= \frac{u_{i+2}^n - 2 u_{i+1}^n +2u_{i-1}^n -u_{i-2}^n}{2\Delta x^3} \end{align*}
To add an interesting fact: this can't be solved via step by step approach. At least as far as i know. BTCS and FTCS methods do not work. Also other schemes are not stable (eigenvalues of step by step matrix form a cross in complex plane and are outside stability regions of the schemes for any $\Delta t$ and $\Delta x$). This is why i would like to use the approach of solving system of equations to get values at grid points.
Here is my try deriving system of equation when i introduce vector $v$. I define $u^{n}:=[u_1^{n},u_2^{n},\dots,u_M^{n}]^T$ and $v^{n}:=[v_1^{n},v_2^{n},\dots,v_M^{n}]^T$ and then $w^n= [(u^n)^T,(v^n)^T]^T$ and then discretizing the space yields
$$ \begin{pmatrix} u^n \\ v^n \end{pmatrix}_t = \begin{pmatrix} 0 & I \\ A & 0 \end{pmatrix}w^n =: D w^n $$ where $$ A=\frac{1}{2 \Delta x^3} \begin{pmatrix} 0 & -2 & 1 & 0 & \cdots & 0 & -1 & 2 \\ 2 & 0 & -2 & 1 & 0 & \cdots & 0 & -1\\ -1 & 2 & 0 & -2 & 1 & 0 & \cdots & 0 \\ 0 & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & 0 \\ 0 & \cdots & 0 & -1 & 2 & 0 & -2 & 1 \\ 1 & 0 & \cdots & 0 & -1 & 2 & 0 & -2 \\ -2 & 1 & 0 & \cdots & 0 & -1 & 2 & 0 \end{pmatrix} $$ using periodic boundary conditions. Central discretization in time yields
$$w^{n+1} - w^{n-1} - 2\Delta tD w^n = 0$$
In the last step i use backward in time derivative approximation to get
$$-w^{N-1} - (\Delta tD - I) w^{N} = 0$$
Lastly, defining $w=[(w^0)^T,(w^1)T,\dots,(w^N)^T]^T$,
\begin{align*} C &=\begin{pmatrix} I & 0 & 0 & 0 & \dots & 0 & 0 & 0 \\ -I & -2\Delta tD & I & 0 & \dots & 0 & 0 & 0 \\ 0 & -I & -2\Delta tD & I & \ddots & 0 & 0 & 0 \\ 0 & 0 & -I & -2\Delta tD & \ddots & \ddots & 0 & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & \ddots & \ddots & \vdots \\ 0 & 0 & 0 & \ddots & \ddots & -2\Delta tD & I & 0 \\ 0 & 0 & 0 & 0 & \ddots & -I & -2\Delta tD & I \\ 0 & 0 & 0 & 0 & \dots & 0 & -I & I-\Delta tD \\ \end{pmatrix}\\ \end{align*} and $c = ((w^0)^T,0...,0)^T$ i have to solve
$$C w = c$$
Edit: It works in Matlab now and the results are nice. Just matrix inversion is not fast.