I'm trying to solve the pde $u_t = i[\begin{pmatrix}0 & 1 \\ 1 & 0\end{pmatrix}u_{xx}+\begin{pmatrix}3 & -1 \\ -1 & 3\end{pmatrix}u]$ where $0 \leq x \leq 1$, $0 \leq t$
with periodic boundary conditions, subject to the initial condition $u(0, x) = \begin{pmatrix}\cos(2\pi x)\\-i\sin(2\pi x) \end{pmatrix}$
Using finite difference forward scheme, also known as forward Euler scheme.
What I tried and my results:
For simplicity suppose $u(t,x) = \begin{pmatrix}f(t,x) \\ g(t,x)\end{pmatrix}$, then our system is $\begin{pmatrix}f_t\\g_t\end{pmatrix} = i\begin{pmatrix}g_{xx}+3f-g\\f_{xx}-f+3g\end{pmatrix}$
Our approximations for the derivatives are:
$f_{t}(t_j,x_i) \approx \frac{f(t_{j+1},x_i)-f(t_j,x_i)}{\Delta t}$, $f_{xx}(t_j,x_i) \approx \frac{f(t_j,x_{i+1})-2f(t_j,x_i)+f(t_j,x_{i-1})}{\Delta x^2}$ and the same for $g$, where $(t_j,x_i)$ are grid points.
Let's look at our first equation, the one for $f_t$. Plugging in the derivative approximations we get $f(t_{j+1},x_i) = f(t_j, x_i )+i(\frac{\Delta t}{\Delta x^2}(g(t_j,x_{i+1})-2g(t_j,x_i)+g(t_j,x_{i-1}))+3\Delta t f-\Delta tg)$
In simpler, shorter matrix form:
$$\begin{pmatrix}f(t_{j+1},0)\\f(t_{j+1},\Delta x)\\ \vdots\\f(t_{j+1},1)\end{pmatrix}=\begin{pmatrix}f(t_{j},0)\\f(t_{j},\Delta x)\\ \vdots\\f(t_{j},1)\end{pmatrix}+i\begin{pmatrix}-\frac{2\Delta t}{\Delta x^2}&\frac{\Delta t}{\Delta x^2} & 0 & \dots & 0 & \frac{\Delta t}{\Delta x^2}\\ \frac{\Delta t}{\Delta x^2}& -\frac{2\Delta t}{\Delta x^2}&\frac{\Delta t}{\Delta x^2}& \dots & 0 & 0\\ \vdots & \vdots & \vdots & \dots & \vdots & \vdots \\\frac{\Delta t}{\Delta x^2} & 0 & 0 & \dots & \frac{\Delta t}{\Delta x^2} & -\frac{2\Delta t}{\Delta x^2}\end{pmatrix}\begin{pmatrix}g(t_j,0) \\ g(t_j, \Delta x) \\ \vdots \\ g(t_j, 1)\end{pmatrix} + 3i\Delta t\begin{pmatrix}f(t_j, 0) \\ f(t_j, \Delta x) \\ \vdots \\ f(t_j, 1)\end{pmatrix}-i\Delta t\begin{pmatrix}g(t_j, 0) \\ g(t_j, \Delta x) \\ \vdots \\ g(t_j, 1)\end{pmatrix}$$
And also the respective system for $g$ at time $t_{j+1}$. And that's exactly what I did. I put in the values we know for $t = 0$, and then i do one step (forward in time) in $f$, one step in $g$, and so on, until $t=1$.
My issue is that my solution seems to explode. The result I get for $f$ at $t=1$ oscillates wildly between $10^{160}$ to $-10^{160}$, and the same for $g$. For this I used $\Delta x = \frac{1}{16}$ and $\Delta t = 0.25 \cdot \Delta x^2$.
Does this result make sense? where is the problem? What is the actual solution?