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. 
Boundary value method for equation $u_{tt} = u_{xxx}$
pde;discretization