## 2 Proposed Numerical Scheme
### Stability Analysis for System of Equations  
The system of linear or linearized first-order differential equations can be expressed in a single vector matrix of the form  
\[\mathbf{v}^{\prime}=A\mathbf{v} \tag{20}\]  
The first or predictor stage of the scheme for the solution of Equation (20) can be expressed as  
\[\overline{v}_{i}=\mathbf{v}_{i-1}+hA\mathbf{v}_{i}+h^{2}A^{2}\mathbf{v}_{i-1} \tag{21}\]  
The second or corrector stage of the scheme for the solution of Equation (20) is given as  
\[\mathbf{v}_{i}=\mathbf{v}_{i-1}+h\left\{\frac{5}{9}A\mathbf{v}_{i}+\frac{1}{2}A\mathbf{v}_{i- 1}-\frac{1}{18}A\overline{\mathbf{v}}_{i}\right\} \tag{22}\]  
Since the proposed scheme is implicit in solving difference Equations (21) and (22), an iterative scheme is employed. Therefore, applying the Gauss-Seidel iterative method to Equation (21) is obtained  
\[\overline{\mathbf{v}}_{i}^{k+1}=\mathbf{v}_{i-1}^{k+1}+hA\mathbf{v}_{i}^{k}+h^{2}A^{2}\mathbf{ v}_{i-1}^{k+1} \tag{23}\]  
Similarly, Equation (22) is re-written as  
\[\mathbf{v}_{i}^{k+1}=\mathbf{v}_{i-1}^{k+1}+h\left\{\frac{5}{9}A\mathbf{v}_{i}^{k}+\frac{1 }{2}A\mathbf{v}_{i-1}^{k+1}-\frac{1}{18}A\overline{\mathbf{v}}_{i}^{k+1}\right\} \tag{24}\]  
The Fourier series or Von Neumann stability analysis is employed to check the proposed scheme's stability conditions. To do so, consider the following transformations:  
\[\overline{\mathbf{v}}_{i}^{k+1}=\mathbf{E}^{k+1}e^{il\psi},\ \mathbf{v}_{i-1}^{k+1}=\mathbf{E}^{k+1}e^{(i-1) l\psi} \tag{25}\]  
\[\mathbf{v}_{i}^{k}=\mathbf{E}^{k}e^{il\psi},\ \mathbf{v}_{i}^{k+1}=\mathbf{E}^{k+1}e^{il\psi}\]  
where \(l=\sqrt{-1}\).  
Substituting some of the relevant transformations from (25) into Equation (21) yields  
\[\overline{\mathbf{E}}^{k+1}e^{il\psi}=\mathbf{E}^{k+1}e^{(i-1)l\psi}+hA\mathbf{E}^{k}e^{ il\psi}+h^{2}A^{2}\mathbf{E}^{k+1}e^{(i-1)l\psi} \tag{26}\]  
Simplifying Equation (26), it is obtained:  
\[\overline{\mathbf{E}}^{k+1}=\mathbf{E}^{k+1}e^{-l\psi}+hA\mathbf{E}^{k}+h^{2}A^{2}\mathbf{E}^ {k+1}e^{-l\psi} \tag{27}\]Similar to the predictor substituting relevant transformation from Equation (25) into Equation (24), the resulting equation is given by  
\[\mathbf{E}^{k+1}e^{i\psi}=\mathbf{E}^{k+1}e^{(i-1)I\psi}+h\bigg{\{}\frac{5}{9}A\mathbf{E}^{k} e^{iI\psi}+\frac{1}{2}A\mathbf{E}^{k+1}e^{(i-1)I\psi}-\frac{1}{18}A\overline{\mathbf{E}}^{k+1 }e^{iI\psi}\bigg{\}} \tag{28}\]  
Upon simplifying Equation (28), it leads to  
\[\mathbf{E}^{k+1}=\mathbf{E}^{k+1}e^{-I\psi}+h\bigg{\{}\frac{5}{9}A\mathbf{E}^{k}+\frac{1}{2} A\mathbf{E}^{k+1}e^{-I\psi}-\frac{1}{18}A\overline{\mathbf{E}}^{k+1}\bigg{\}} \tag{29}\]  
Inserting Equation (27) into Equation (29) gives  
\[\mathbf{E}^{k+1}=\mathbf{E}^{k+1}e^{-I\psi}+h\bigg{\{}\frac{5}{9}A\mathbf{E}^{k}+\frac{1}{2 }A\mathbf{E}^{k+1}e^{-I\psi}-\]  
Re-write Equation (30) as  
\[\bigg{[}I.D-I.De^{-I\psi}-\frac{h}{2}Ae^{-I\psi}+\frac{h}{18}Ae^{-I\psi}+\frac {h^{3}}{18}A^{3}e^{-I\psi}\bigg{]}\mathbf{E}^{k+1}=\bigg{(}\frac{5}{9}hA-\frac{h^ {2}A^{2}}{18}\bigg{)}\mathbf{E}^{k} \tag{31}\]  
The stability can be expressed as  
\[\Bigg{|}\frac{\frac{5h\lambda_{A}}{9}-\frac{h^{2}}{18}\lambda_{A}^{2}}{1-e^{ -I\psi}-\frac{h}{2}\lambda_{A}e^{-I\psi}+\frac{h}{18}\lambda_{A}e^{-I\psi}+ \frac{h^{3}\lambda_{A}^{3}}{18}e^{-I\psi}}\Bigg{|}<1 \tag{32}\]  
where \(I.D.\) is the identity matrix.