I'd like to solve the heat equation numerically: \begin{equation} \partial_tu=\partial_x^2u, \end{equation} and I've tested my algorithm with initial conditions $u(0,x)\propto\operatorname{exp}[-(x-x_0)/(2\sigma^2)]$ and boundary conditions $u(t,0)=u(t,1)=0$. However I'm not obtaining stable regimes as I would expect from von-Neumann stability analysis. Where could this instability come from?
Problem description: Let $u^n_j=u(n\Delta t,j\Delta x)$ denote $u$ on a grid. When discretising the time derivative according to $\partial_tu=\frac{1}{2\Delta t}(u^{n+1}_j-u^{n-1}_j)$, I yield the result that I'll always have unstable modes. This is what I'm actually observing.
If I instead apply the Euler method (an actually less exact approximation) to the time derivative, $\partial_tu=\frac{1}{\Delta t}(u^{n+1}_j-u^n_j)$, I yield the stable regime \begin{equation} \Delta x^2\geq\Delta t. \end{equation} However, if I exactly fulfil the equality $\Delta x^2=\Delta t$ (e. g., $\Delta x=10^{-3},\Delta t=10^{-6}$), I do not yield a stable result. I achieve longer stability by further increasing $\Delta x^2>\Delta t$ and for $\Delta x\geq1.5\sqrt{\Delta t}$ I haven't observed instability so far.
What am I doing wrong? Where does the instability come from? According to the stability analysis, I shouldn't have any exponentially growing modes at all as long as I choose the step sizes appropriately.
My thoughts: There are three sources of error: Rounding, discretisation of the derivatives, subtractive cancellation. Rounding errors are of arbitrary wave length, but small in amplitude. Since no mode is amplified (stability criterion), these errors shouldn't lead to divergence of my solution. Could any of the other two error sources be the cause of divergence?