Let's pretend we have a spatially discretized PDE of the following form:\begin{align}\frac{\partial^2 \boldsymbol{u}^k}{\partial t^2} = D\boldsymbol{u}^k\end{align}where $D$ can be any form for now and $k$ refers to this being the discretization at time $t_k$. Then let's suppose, for example, we use a central difference approximation for the time derivative. We would then arrive at the following:\begin{align}\frac{\partial^2 \boldsymbol{u}^k}{\partial t^2} &= D\boldsymbol{u}^k \\\frac{\boldsymbol{u}^{k+1} - 2\boldsymbol{u}^k + \boldsymbol{u}^{k-1}}{\Delta t^2} &= D\boldsymbol{u}^k \\\boldsymbol{u}^{k+1} &= \left( 2I + \Delta t^2 D \right)\boldsymbol{u}^k - \boldsymbol{u}^{k-1}\\\end{align}Given the formulation above, how could one approach stability? First thoughts go to using Von Neumann stability analysis, but I have never seen it used in a vector-wise fashion such that one could take into account an arbitrary $D$ matrix. Any thoughts or references would be very useful.EditThe references provided in the comments were useful, but I found that only a few simple connections were needed to approach stability for this scenario. The key link was to take the expression above and cast it into a different difference equation with a more suitable form. The more suitable case is by defining $\boldsymbol{w}^k = \left[(\boldsymbol{u}^{k})^{T}, (\boldsymbol{u}^{k-1})^{T}\right]^T$. We can than recast our difference equation into the form:\begin{align}\boldsymbol{w}^{k+1} = G \boldsymbol{w}^k\end{align}where \begin{align}G = \begin{bmatrix}(2I + \Delta t^2 D) & -I \\I & 0\end{bmatrix}\end{align}Then we know this system is stable if, given the set of eigenvalues $\lbrace \lambda_i \rbrace$ for $G$, the following holds:\begin{align}\max_{i} \left| \lambda_i\right| \lt 1\end{align}With these results, we can check whether some differential operator $D$ is going to be stable with a given time discretization, particularly if it uses older $\boldsymbol{u}$ states.
Stability of PDE Discretizations with Multistep Time Discretizations
pde;numerical analysis;stability;differential equations
What you're looking for is a discussion of absolute stability of multistep methods for ODEs. See for example chapter 7 of LeVeque's book on finite difference methods . Once you know the region of absolute stability for your time discretization, it's just a matter of checking that the eigenvalues of D (scaled by the time step) are inside that region. This is typical method-of-lines stability analysis. The approach you've mentioned in your edit is equivalent, but in a sense less flexible since it doesn't separate the analysis of the time and space discretizations so clearly.   Also, the approach I'm recommending can be done entirely with pencil and paper. 