## 2 Results
### Numerical Stability  
It is readily verified that the the solution to (25) is  
\[w(\tau_{n+1})=e^{\Delta\tau A}w(\tau_{n})+\int_{\tau_{n}}^{\tau_{n+1}}e^{(t-\tau_ {n})A}f(t)dt,\;\;n=0,1,\ldots, \tag{29}\]  
where \(\tau_{n}=n\Delta\tau\). The formal solution to (28) is similar. We have:  
**Lemma 1** ([10; 19]). : _The semi-discretized schemes (25) and (28) are stable if_  
\[\lim_{h_{max},k_{max}\to 0}\left(\max_{\tau\in[0,\tau^{*}]}\left\|e^{\tau A} \right\|_{2}\right)\leq c(\tau^{*}),\;\lim_{h_{max},k_{max}\to 0}\left(\max_{ \tau\in[0,\tau^{*}]}\left\|e^{\tau A}\right\|_{2}\right)\leq c(\tau^{*}),\]  
_where \(\tau^{*}\in(0,T)\)._  
**Lemma 2** ([10]). : _Let \(B\in\mathbb{C}^{d\times d}\). Then \(\sigma(B)\subset\cup_{i=1}^{d}S_{i}\), where_  
\[S_{i}=\left\{z\in C:|z-b_{i,j}|\leq\sum_{j=1,j\
eq i}^{d}|b_{i,j}|\right\}\]  
_are Gershgorin discs, and \(\sigma(B)\) is the set of all eigenvalues of \(B\). Moreover, \(\lambda\in\sigma(B)\) may lie on \(\partial S_{\partial}\) for some \(i^{0}\in\{1,2,\ldots,d\}\), only if \(\lambda\in\partial S_{i}\) for all \(i=1,2,\ldots,d\)._  
**Lemma 3** ([20]). : _The matrix exponential \(e^{tA}\) tends to a zero matrix as \(t\to+\infty\) if and only if all the eigenvalues of \(A\) have strictly negative real parts._  
**Theorem 1**. : _The semi-discretized schemes (25) and (28) are linearly stable._  
Proof. : We will only need to show the case of \(\rho\in(0,1],\;\eta>2r\) for (25), since extensions of our results for other cases are technically imminent. Thus, we only need to show that each of the \(MN\) Gershgorin discs of \(A\) lies on the left side of the complex plane. In fact, there are five types of the Gershgorin discs to consider:  
1. Discs centered at an internal mesh point;
2. discs centered on one of the Dirichlet boundaries;
3. discs centered on the Neumann boundary;
4. discs centered at one of the intersection mesh points of two Dirichlet boundaries; and
5. discs centered at one of the intersection mesh points of one Dirichlet boundary and the Neumann boundary. We provide detailed proofs for the first three types of discs. Similar arguments can be applied to the rest of the cases. Case 1: In this situation, we first consider the situation in which \(\eta<y_{n}\leq Y\). Let \(z\in S_{i}\) be any complex number, where \(S_{i}\) is a Gershgorin disc centered at an internal point of the spatial grids. Thus,  
\[ \begin{split}&\left|z+\frac{y_{n}}{h_{m}h_{m+1}}+\frac{\sigma^{2}y _{n}}{k_{n}k_{n+1}}-\frac{\rho\sigma y_{n}}{2h_{m+1}k_{n+1}}-\frac{\rho \sigma y_{n}}{2h_{m}k_{n}}+\frac{y_{n}-2r}{2h_{m}}-\frac{\kappa(\eta-y_{n})}{ k_{n}}\right|\\ &\leq\left|\frac{\sigma^{2}y_{n}}{k_{n}(k_{n}+k_{n+1})}-\frac{ \rho\sigma y_{n}}{2h_{m}k_{n}}-\frac{\kappa(\eta-y_{n})}{k_{n}}\right|+\left| \frac{\sigma^{2}y_{n}}{k_{n+1}(k_{n}+k_{n+1})}-\frac{\rho\sigma y_{n}}{2h_{m +1}k_{n+1}}\right|\\ &\qquad\qquad+\left|\frac{y_{n}}{h_{m+1}(h_{m}+h_{m+1})}-\frac{ \rho\sigma y_{n}}{2h_{m+1}k_{n+1}}\right|+\left|\frac{\rho\sigma y_{n}}{2h_{m +1}k_{n+1}}\right|+\left|\frac{\rho\sigma y_{n}}{2h_{m}k_{n}}\right|\\ &\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad \qquad+\left|\frac{y_{n}}{h_{m}(h_{m}+h_{m+1})}-\frac{\rho\sigma y_{n}}{2h_{m }k_{n}}+\frac{y_{n}-2r}{2h_{m}}\right|.\end{split}  \tag{30}\]Let \(\alpha\) be the real part of \(z\). Since we are concerned only about the upper bound of the real part of the eigenvalues, we may replace \(z\) by \(\alpha\) via a triangle inequality, and remove the absolute value sign on the left hand side of (30).