I am using Heun's method with a third order upwind spatial scheme, which is suggested by Shao (2008) to be used for solving the horizontal advection part of the advection-diffusion equation.
This is what I got:
\begin{align} &C^{\ast} = C^{n} - A^{n} \delta t\\ &C^{n+1} = C^{n} - \dfrac{1}{2}\left( A^{n} + A^{\ast} \right) \delta t \end{align}
Assuming $u=cte>0$, we have,
\begin{equation} A = \dfrac{u}{\delta x} ( \dfrac{1}{6} C_{i-2} - C_{i-1} + \dfrac{1}{2} C_{i} + \dfrac{1}{3} C_{i+1} ) \end{equation}
\begin{equation} \hat{C}_{\ast} = \hat{C}_{n} - \dfrac{u \delta t}{\delta x} ( \dfrac{1}{6} e^{-2imh} - e^{-imh} + \dfrac{1}{2} + \dfrac{1}{3} e^{imh} ) \hat{C}_{n} \end{equation}
\begin{equation} \hat{C}_{\ast} = \left[ 1 - r F(h) \right] \hat{C}_{n} \end{equation}
In which,
\begin{equation} r = \dfrac{u \delta t}{\delta x} \end{equation} \begin{equation} F(h) = \dfrac{1}{2}+\dfrac{1}{6} e^{-2imh} - e^{-imh} + \dfrac{1}{3} e^{imh} \end{equation}
\begin{equation} \hat{C}_{n+1} = \hat{C}_{n} - \dfrac{r}{2} F(h) \hat{C}_{n} - \dfrac{r}{2} F(h) \hat{C}_{\ast} \end{equation}
\begin{equation} g = \dfrac{\hat{C}_{n+1}}{\hat{C}_{n}} = 1 - r F(h) - \dfrac{r^2}{2} F^2(h) \end{equation}
Enforcing \begin{equation} |g| \leq 1 \end{equation}
should give the stability conditions. I have drawn $g$'s graph
and it seems the scheme is always unstable. But I know it is not. What am I doing wrong? Do you have any suggestions?
Shao, Yaping. Physics and modelling of wind erosion. Springer, 2008.