I'm working on solving the coupled one-dimensional poroelasticity equations (biot's model), given as:$$-(\lambda+ 2\mu) \frac{\partial^2 u}{\partial x^2} + \frac{\partial p}{\partial x} = 0$$$$\frac{\partial}{\partial t} \left[ \gamma p + \frac{\partial u}{\partial x}\right] -\frac{\kappa}{\eta}\left[\frac{\partial^2 p}{\partial x^2}\right]  =q(x,t)$$on the domain $\Omega=(0,1)$ and with the boundary conditions:  $p=0, (\lambda + 2\mu)\frac{\partial u}{\partial x}=-u_0$ at $x=0$ and $u=0, \frac{\partial p}{\partial x} = 0$ at $x=1$.I discretized these equations using a centered finite difference scheme:  $$(\lambda + 2\mu)\frac{u^{t+1}_{i+1}-2u^{t+1}_i+u^{t+1}_{i-1}}{\Delta x^2} + \frac{p^{t+1}_{i+1}-p^{t+1}_{i-1}}{2\Delta x} = 0$$$$\gamma \frac{p^{t+1}_i-p^t_i}{\Delta t} + \frac{u^{t+1}_{i+1}-u^{t+1}_{i-1}}{2\Delta x \Delta t} -\left[ \frac{u^t_{i+1}-u^t_{i-1}}{2\Delta x \Delta t}\right] - \frac{\kappa}{\eta}\left[ \frac{p^{t+1}_{i+1} -2p^{t+1}_i + p^{t+1}_{i-1}}{\Delta x^2}\right]=q^{t+1}_i$$I'm currently working out the details of the scheme's convergence by analyzing its consistency and stability.  The consistency part seems fairly straightforward to me, but I'm already foreseeing some difficulties with the stability analysis.  First of all, there are two variables and two equations.  Secondly, there is also a mixed spatiotemporal derivative term in the second equation.  I'm familiar with von neumann stability analysis and can see that it will be very tough to establish stability with this method.  Are there any alternatives to von neumann analysis that I can use?
Alternatives to von neumann stability analysis for finite difference methods
pde;finite difference;stability;numerical analysis