Non-polynomial spline approach for solving system of singularly perturbed delay differential equations of large delay

ABSTRACT

A computational technique for analysing coupled system of singularly perturbed convection-diffusion delay differential equations with large delay is presented in this study. We solve this problem using non-polynomial spline technique and arrive at a system that can potentially be solved. To show the theoretical investigation, an example is solved for different perturbation parameter values, and the computational results are shown in tables. The model we have proposed is of first order convergent, and this technique indicated that the acquired findings were more accurate than earlier results.

KEYWORDS:

• Singuarly perturbed delay differential equation
• non-polynomial spline
• system of convection diffusion equations

1. Introduction

A singularly perturbed delay differential equation is one in which the highest order derivative term is multiplied by a positive small parameter and at least one delay or advance term, or both, is involved. The study of singular perturbations has come a long way in mathematics and has important potential applications in many areas of science and technology. Different types of real-world phenomena such as in the modelling of the pupil-light reflex in humans (Longtin and Milton 1988), predator-prey model (Gourley and Kuang 2004), HIV infection model (Culshaw and Ruan 2000), variational problem in control theory (Glizer 2000), determination of the behaviour of a neuron to random synaptic input (Lange and Miura 1982), optically bistable device (Mayer et al. 1995),the dynamics of two identical amplifier networks (Chen and Wu 2001), the first exit-time problem (Kot 2001), hydrodynamics of liquid helium (Kuang 1993), various physiological process models (Mackey and Glass 1977) and numerous other applied mathematical problems can be modelled mathematically by singularly perturbed delay differential equations. These differential equations are used to solve problems in which the future is dependent on the present as well as past states of the system being investigated under consideration.

In singular perturbation problems, the convection-diffusion type is referred to when the perturbation parameter equals zero and the order is lowered by one, whereas the reaction-diffusion type is referred to when the order is lowered by two. As a result, second-order singularly perturbed delay differential equations can be either convection-diffusion or reaction-diffusion in nature. In the first type of problem, the left or right-boundary layer is only determined by the sign of the coefficients in the diffusion and convection terms. On the other hand, in reaction-diffusion delay differential equations, there are two boundary layers. As a result, numerous researchers have been working to create various numerical techniques for handling these kinds of problems. When normal numerical methods are used, it is widely known that these types of problems yield inaccurate results for small values of the perturbation parameter. As a result, it is essential to develop numerical techniques that can provide improved accuracy regardless of parameter values. Another crucial aspect of these types of problems is the investigation of the effect of shift parameters on the solution’s boundary layer behaviour. Consequently, there has been a rise in interest in numerical techniques for singularly perturbed delay differential equations in recent years.

In problems involving control systems, the state of the system is observed and adjustments are incorporated to control the outcomes. Due to the impossibility of instantaneous adjustment, there is always a delay between observation and control action. Differential equations with delay arguments govern this type of system. For a single singularly perturbed delay differential equation, robust numerical techniques have been developed by Lange and Miura (Lange and Miura 1994) investigated a class of boundary-value problems and discussed an asymptotic method to approximately solve this type of differential equations. Phaneendra and Lalu (Phaneendra and Lalu 2019) had carried out a quadrature technique for the solution of this type of problems. They used a neutral delay differential equation which is equivalent to the SPDDE. Chakravarthy and Kumar (Chakravarthy and Kumar 2021) proposed Numerov’s method which uses a fitted operator finite difference technique. Sekar and Tamilselvan (Erdogan et al. 2020) suggested a piecewise Shishkin type mesh to solve the problem by using finite difference method. Erdogan and Sakar (Woldaregay and Liu 2022) suggested a numerical technique based on a piecewise uniform Shishkin mesh with an exponentially fitted difference scheme for each time-subinterval. Woldaregay and Durressa (Sekar and Tamilselvan 2019) used Taylor series approximation to approximate the delay terms, and the obtained singularly perturbed boundary value problem is solved with an exponentially fitted operator mid-point upwind finite difference method. But for systems of equations, only a small number of results have been reported. Subburayan and Ramanujam (Subburayan and Ramanujam 2014) proposed an asymptotic computational technique for convection diffusion type equations. The authors in (Selvi and Ramanujam 2017) applied a computational technique to solve a system of singularly perturbed delay differential equations of large delay that was weakly coupled. Chakravarthy and Gupta (Chakravarthy and Gupta 2020) focus on the use of cubic splines in tension approximation techniques to solve system of singularly perturbed convection diffusion equation.

We developed a fitted non-polynomial spline technique for the numerical solution of the system of second-order singularly perturbed delay differential equations with large delay. When the perturbation parameter is small compared to the mesh width used to discretize the given problem, it is well known that the conventional techniques fail. Our goal is to demonstrate that, when epsilon is small or large in relation to step size h, accurate numerical approximations can be obtained using a non-polynomial spline. Here, we take a non-polynomial spline approach using the fitted operator method to solve this kind of problem. In Sect. 2, we define the problem and assumptions on the parameter $\mathit{\epsilon }$ and stability analysis is discussed in Sect. 3. In Sect. 4 we employ the non-polynomial spline to arrive at its numerical solution by Gauss elimination method. In Sect. 5, we perform the error and convergence analysis of the method. In Sect. 6 we took some examples and compared the results with (Subburayan and Ramanujam 2014; Selvi and Ramanujam 2017; Chakravarthy and Gupta 2020) and their results are shown. Conclusion follows in Sect. 7.

2. Statement of the problem

Consider the following class of system of singularly perturbed delay differential equation on the domain $\mathrm{\Gamma }=(0,2)$,

(1) $\left.\begin{array}{c}{\mathcal{L}}_1\mathbf{z}\mathbf{(}\mathbf{t}\mathbf{)}\equiv -\mathit{\epsilon }{\mathit{z}{\mathit{ }}^{\mathrm{\prime \prime }}}_1(\mathit{t})+p_1(\mathit{t}){\mathit{z}{\mathit{ }}^{\mathrm{\prime }}}_1(\mathit{t})+{\sum }_{k=1}^2{q}_{1\mathit{k}}(\mathit{t})z_k(\mathit{t})+{\sum }_{k=1}^2{r}_{1\mathit{k}}(\mathit{t})z_k(\mathit{t}-1)=s_1(\mathit{t}),\mathit{t}\in (0,2)\\ {\mathcal{L}}_2\mathbf{z}\mathbf{(}\mathbf{t}\mathbf{)}\equiv -\mathit{\epsilon }{\mathit{z}{\mathit{ }}^{\mathrm{\prime \prime }}}_2(\mathit{t})+p_2(\mathit{t}){\mathit{z}{\mathit{ }}^{\mathrm{\prime }}}_2(\mathit{t})+{\sum }_{k=1}^2{q}_{2\mathit{k}}(\mathit{t})z_k(\mathit{t})+{\sum }_{k=1}^2{r}_{2\mathit{k}}(\mathit{t})z_k(\mathit{t}-1)=s_2(\mathit{t}),\mathit{t}\in (0,2)\end{array}\right\}$(1)

with boundary conditions,

(2) $\begin{array}{cc}{z}_1(\mathit{t})={\varphi }_1(\mathit{t}),& \mathit{t}\in [-1,0],z_1(2)={\gamma }_1,\\ z_2(\mathit{t})={\varphi }_2(\mathit{t}),& \mathit{t}\in [-1,0],z_2(2)={\gamma }_2,\end{array}$(2)

where the perturbation parameter satisfy $0<\mathit{\epsilon }<<1$.

It is also assumed that the functions $p_1(\mathit{t}),p_2(\mathit{t}),q_{11}(\mathit{t}),$$q_{12}(\mathit{t}),q_{21}(\mathit{t}),q_{22}(\mathit{t}),\text{break}{r}_{11}(\mathit{t}),r_{12}(\mathit{t}),r_{21}(\mathit{t}),r_{22}(\mathit{t}),s_1(\mathit{t})$ and $s_2(\mathit{t})$ are sufficiently smooth functions on $\bar{\mathrm{\Gamma }}=[0,2]$ together with the following conditions:

$p_1(\mathit{t})\ge l_1>0\mathrm{and}{p}_2(\mathit{t})\ge l_2>0$

and

$\underset{t\in \bar{\mathrm{\Gamma }}}{min}\left\{\sum _{\mathit{k}=1}^2{q}_{1\mathit{k}}(\mathit{t}),\sum _{\mathit{k}=1}^2{q}_{2\mathit{k}}(\mathit{t})\right\}\ge {\lambda }_1>0,$

$\underset{t\in \bar{\mathrm{\Gamma }}}{min}\left\{\sum _{\mathit{k}=1}^2{r}_{1\mathit{k}}(\mathit{t}),\sum _{\mathit{k}=1}^2{r}_{2\mathit{k}}(\mathit{t})\right\}\ge {\lambda }_2>0.$

The differential operator in (1) can be written as:

$\mathcal{L}\mathbf{z}={\left({\mathcal{L}}_1\mathbf{z},{\mathcal{L}}_2\mathbf{z}\right)}^{\mathit{T}}=-\mathit{\epsilon }\frac{d^2\mathbf{z}}{\mathit{d}{t}^2}+\mathit{P}(\mathit{t})\frac{\mathit{d}\mathbf{z}}{\mathit{dt}}+\mathit{Q}(\mathit{t})\mathbf{z}+\mathit{R}(\mathit{t})\mathbf{z}(\mathit{t}-1),$

where the coefficient matrices are given by,

$\mathit{\epsilon }=\left[\begin{array}{cc}\mathit{\epsilon }& 0\\ 0& \mathit{\epsilon }\end{array}\right],\mathit{P}(\mathit{t})=\left[\begin{array}{cc}{p}_1(\mathit{t})& 0\\ 0& p_2(\mathit{t})\end{array}\right],\mathit{Q}(\mathit{t})=\left[\begin{array}{cc}{q}_{11}& q_{12}\\ q_{21}& q_{22}\end{array}\right],\mathit{R}(\mathit{t})=\left[\begin{array}{cc}{r}_{11}& r_{12}\\ r_{21}& r_{22}\end{array}\right].$

Also the functions $p_i,q_{\mathit{ik}},r_{\mathit{ik}},s_i\in C^4(\bar{\mathrm{\Gamma }}),\mathit{i}=1,2,\mathit{k}=1,2$. In addition, the boundary conditions $\mathit{\varphi }={\left({\varphi }_1,{\varphi }_2\right)}^{\mathit{T}},\mathit{\gamma }={\left({\gamma }_1,{\gamma }_2\right)}^{\mathit{T}}$ and the source term $\mathbf{s}={\left(s_1,s_2\right)}^{\mathit{T}}$ are smooth functions on ${\mathrm{\Gamma }}^{\ast }$, where $\mathrm{\Gamma }=(0,2),\bar{\mathrm{\Gamma }}=[0,2],{\mathrm{\Gamma }}^{\ast }=[-1,0]$. In this context, $C^n(\mathrm{\Gamma })$ refers to the class of $\mathit{n}$ times continuously differentiable functions in $\mathrm{\Gamma }$. Also, the problem (1) and (2) shows a weak boundary layer at $\mathit{t}=1$ and a strong layer at $\mathit{t}=2$.

3. Stability analysis

Theorem 3.1.

Let the function $\mathbf{z}={\left(z_1,z_2\right)}^{\mathit{T}},z_1,z_2\in C^0\left(\bar{\mathrm{\Gamma }}\right)\cap C^2({\mathrm{\Gamma }}^{\ast })$ satisfying $z_j(0)\ge 0,z_j(2)\ge 0,{\mathcal{L}}_j\mathbf{z}\ge 0,\mathit{j}=1,2,\mathrm{\forall t}\in {\mathrm{\Gamma }}^{\ast }$ and $\mathit{z}{{\mathit{ }}^{\mathrm{\prime }}}_{\mathit{j}}(1+)-z_{\mathit{j}{\mathit{ }}^{\mathrm{\prime }}}(1-)=[\mathit{z}{{\mathit{ }}^{\mathrm{\prime }}}_{\mathit{j}}](1)\le 0,\mathit{j}\text{break}=1,2$. Then $z_j(\mathit{t})\ge 0,\mathrm{\forall t}\in \bar{\mathrm{\Gamma }},\mathit{j}=1,2.$

Proof.

Define $\mathbf{r}=(r_1,r_2)$, where

$r_1(\mathit{t})=r_2(\mathit{t})=\left\{\begin{array}{ccc}\frac{1}{8}+\frac{t}{2}& ,& \mathit{t}\in [0,1]\\ \frac{3}{8}+\frac{t}{4}& ,& \mathit{t}\in [1,2]\end{array}\right.$

Here $r_i(\mathit{t})$ is positive for all $\mathit{x}\in \bar{\mathrm{\Gamma }},{\mathcal{L}}_i\mathbf{r}>0,\mathrm{\forall }\mathit{t}\in {\mathrm{\Gamma }}^{\ast }$ and $[\mathit{r}{{\mathit{ }}^{\mathrm{\prime }}}_{\mathit{i}}](1)<0,\mathit{i}=1,2$. Let

$\mathit{\mu }=max\left\{\underset{t\in \bar{\mathrm{\Gamma }}}{max}\left\{-\frac{z_1(\mathit{t})}{r_1(\mathit{t})}\right\},\underset{t\in \bar{\mathrm{\Gamma }}}{max}\left\{-\frac{z_2(\mathit{t})}{r_2(\mathit{t})}\right\}\right\}.$

Then there exists $t_0\in \bar{\mathrm{\Gamma }}$ satisfying $z_1(t_0)+\mathit{\mu }{r}_1(t_0)=0$ or $z_2(t_0)+\mathit{\mu }{r}_2(t_0)=0$ or both and $z_j(\mathit{t})+\mathit{\mu }{r}_j(\mathit{t})\ge 0,\mathrm{\forall t}\in \bar{\mathrm{\Gamma }},\mathit{j}=1,2$. Without loss of generality, we assume that $z_1(t_0)+\mathit{\mu }{r}_1(t_0)=0$. Hence $z_1+\mathit{\mu }{r}_1$ attains its minimum. Suppose the theorem does not hold true, then $\mathit{\mu }>0$.

Case 1: $t_0\in {\mathrm{\Gamma }}^{-}=(0,1).$

Then, for $\mathit{j}=1,2$,

$\begin{array}{rl}& 0<{\mathcal{L}}_j(\mathbf{z}+\mathit{\mu }\mathbf{r})(t_0)=-\mathit{\epsilon }(z_j+\mathit{\mu }{r}_j){ }^{\mathrm{\prime \prime }}(t_0)+p_j(t_0)(z_j+\mathit{\mu }{r}_j){ }^{\mathrm{\prime }}(t_0)\\ & +\sum _{\mathit{i}=1}^2{q}_{\mathit{ji}}(t_0)(z_j+\mathit{\mu }{r}_j)(t_0)\le 0,\end{array}$

which is a contradiction to our assumption.

Case 2: $t_0\in {\mathrm{\Gamma }}^{+}=(1,2)$.

Then, for $\mathit{j}=1,2$,

$\begin{array}{rl}& 0<{\mathcal{L}}_j(\mathbf{z}+\mathit{\mu }\mathbf{r})(t_0)=-\mathit{\epsilon }(z_j+\mathit{\mu }{r}_j){ }^{\mathrm{\prime \prime }}(t_0)+p_j(t_0)(z_j+\mathit{\mu }{r}_j){ }^{\mathrm{\prime }}(t_0)\\ & +\sum _{\mathit{i}=1}^2{q}_{\mathit{ji}}(t_0)(z_j+\mathit{\mu }{r}_j)(t_0)+\sum _{\mathit{i}=1}^2{r}_{\mathit{ji}}(t_0)(z_j+\mathit{\mu }{r}_j)(t_0-1)\le 0,\end{array}$

which is a contradiction to our assumption.

Case 3: $t_0=1.$

Then, for $\mathit{j}=1,2$,

$0\le [(z_j+\mathit{\mu }{r}_j){ }^{\mathrm{\prime }}](1)=[\mathit{z}{{\mathit{ }}^{\mathrm{\prime }}}_{\mathit{j}}](1)+\mathit{\mu }[\mathit{r}{{\mathit{ }}^{\mathrm{\prime }}}_{\mathit{j}}](1)<0,$

which is a contradiction to our assumption.

From this, it is clear that in all the three cases we arrived a contradiction. Therefore $\mathit{\mu }>0$ is not possible. This shows that,

$z_j(\mathit{t})\ge 0,\mathrm{\forall t}\in \bar{\mathrm{\Gamma }},\mathit{j}=1,2.$

Theorem 3.2.

Let $\mathbf{z}=(z_1,z_2),z_1,z_2\in C^0(\bar{\mathrm{\Gamma }})\cap C^1(\mathrm{\Gamma })\cap C^2({\mathrm{\Gamma }}^{\ast })$ be any function. Then

$|z_j(\mathit{t})|\le \mathit{K}max\left\{\underset{\mathit{i}=1,2}{max}\left\{|z_i(0)|\right\},\underset{\mathit{i}=1,2}{max}\left\{|z_i(2)|\right\},\underset{\mathit{i}=1,2}{max}\left\{\parallel {\mathcal{L}}_i\mathbf{z}{\parallel }_{{\mathrm{\Gamma }}^{\ast }}\right\}\right\}$

for all $\mathit{t}\in \bar{\mathrm{\Gamma }},\mathit{j}=1,2.$

Proof. Let $\mathit{K}>0$ be a constant. Define ${\tau }^{\pm }=({\tau }_1^{\pm },{\tau }_2^{\pm })$, where,

${\tau }_j^{\pm }(\mathit{t})=\mathit{K}\bar{K}{r}_j(\mathit{t})\pm z_j(\mathit{t}),\mathrm{\forall t}\in \bar{\mathrm{\Gamma }},\mathit{j}=1,2,$

where,

$\bar{K}=max\left\{\underset{\mathit{i}=1,2}{max}\left\{|z_i(0)|\right\},\underset{\mathit{i}=1,2}{max}\left\{|z_i(2)|\right\},\underset{\mathit{i}=1,2}{max}\left\{\parallel {\mathcal{L}}_i\mathbf{z}{\parallel }_{{\mathrm{\Gamma }}^{\ast }}\right\}\right\}.$

By appropriate choice of $\mathit{K}$, we obtain,

${\tau }_j^{\pm }(0)=\mathit{K}\bar{K}{r}_j(0)\pm z_j(0)>0\mathrm{and}{\tau }_j^{\pm }(1)=\mathit{K}\bar{K}{r}_j(1)\pm z_j(1)>0,\mathit{j}=1,2,$

Let $\mathit{t}\in {\mathrm{\Gamma }}^{-}$, then by appropriate choice of $\mathit{K}$ and $\mathit{j}=1,2$, we obtain,

${\mathcal{L}}_j{\tau }^{\pm }=\mathit{K}\bar{K}{\mathcal{L}}_j\mathbf{r}\pm {\mathcal{L}}_j\mathbf{z}\ge 0.$

Similarly we can prove for ${\mathcal{L}}_j{\tau }^{\pm }\ge 0$ in ${\mathrm{\Gamma }}^{+}$. Therefore ${\mathcal{L}}_j{\tau }^{\pm }\ge 0$ in ${\mathrm{\Gamma }}^{\ast }$. Also by proper choice of $\mathit{K},$ we obtain,

$\left[{\tau }_j^{\pm {\mathit{ }}^{\mathrm{\prime }}}\right](1)=\mathit{K}\bar{K}[\mathit{r}{{\mathit{ }}^{\mathrm{\prime }}}_{\mathit{j}}](1)+[\mathit{z}{{\mathit{ }}^{\mathrm{\prime }}}_{\mathit{j}}](1)<0,\mathit{j}=1,2.$

Then by Theorem 3.1, we obtain ${\tau }_j^{\pm }(\mathit{t})\ge 0,\mathit{t}\in \bar{\mathrm{\Gamma }},\mathit{j}=1,2$. Therefore,

$|z_j(\mathit{t})|\le \mathit{K}max\left\{\underset{\mathit{i}=1,2}{max}\left\{|z_i(0)|\right\},\underset{\mathit{i}=1,2}{max}\left\{|z_i(2)|\right\},\underset{\mathit{i}=1,2}{max}\left\{\parallel {\mathcal{L}}_i\mathbf{z}{\parallel }_{{\mathrm{\Gamma }}^{\ast }}\right\}\right\},$

for all $\mathit{t}\in \bar{\mathrm{\Gamma }},\mathit{j}=1,2.$

4. Derivation of the method

The objective of this section is to introduce the uniform mesh and the numerical approach that is being examined in this article. Let ${\bar{\mathrm{\Gamma }}}^{2\mathit{N}}=\{t_i{\}}_{i=0}^{2N}$ be the discretized domain with $2\mathit{N}$ mesh – intervals on the domain $\bar{\mathrm{\Gamma }}=[0,2]$, where $\mathit{N}$ is as an even positive integer. Let $0=t_0,t_1,\dots ,t_N=1,t_{\mathit{N}+1},\dots ,t_{2\mathit{N}}=2$ be the mesh points obtained while dividing $[0,2]$, such that the step size $h_i=t_i-t_{\mathit{i}-1}$ for $\mathit{i}=1,2,\dots ,2\mathit{N}-1.$

4.1. Numerical algorithm

The following algorithm is proposed to obtain the numerical solution:

Step 1: Introduce the uniform mesh by discretizing the domain $[0,2]$ with $2\mathit{N}$ mesh intervals.

Step 2: Non-polynomial spline method is applied to the statement problem to obtain the scheme.

Step 3: We use Taylor series approximations for the derivatives in the scheme obtained in Step 2 to arrive at two systems for $Z_1$ and $Z_2$.

Step 4: Introduce fitting factor to the schemes obtained in Step 3, and we derive the schemes into a system of equations in different mesh intervals, incorporating history functions.

Step 5: Finally, we employ Gauss elimination method to solve the system obtained in Step 4 by using MATLAB R2022a mathematical software.

4.2. The proposed numerical scheme

The non-polynomial spline $M_{{\mathrm{\Delta }}_j}(\mathit{t})$ is of the following form:

(3) $\begin{array}{rl}& M_{{\mathrm{\Delta }}_j}(\mathit{t})=a_{\mathit{j},\mathit{i}}sin\mathit{\omega }(\mathit{t}-t_i)+b_{\mathit{j},\mathit{i}}cos\mathit{\omega }(\mathit{t}-t_i)+c_{\mathit{j},\mathit{i}}\left(e^{\omega (t-t_i)}-e^{-\omega (t-t_i)}\right)\\ & +d_{\mathit{j},\mathit{i}}\left(e^{\omega (t-t_i)}+e^{-\omega (t-t_i)}\right),\end{array}$(3)

for $\mathit{j}=1,2\mathrm{and}\mathit{i}=1,2,\dots ,2\mathit{N}-1,$ where the coefficients $a_{\mathit{j},\mathit{i}},b_{\mathit{j},\mathit{i}},c_{\mathit{j},\mathit{i}},d_{\mathit{j},\mathit{i}}$ are unknown and $\mathit{\omega }\ne 0$ is an arbitrary parameter that will be utilized to improve the methods accuracy.

In order to calculate the coefficients which are unknown in (3), we use the conditions,

$M_{{\mathrm{\Delta }}_j}(t_i)=Z_{\mathit{j},\mathit{i}},M_{{\mathrm{\Delta }}_j}(t_{\mathit{i}+1})=Z_{\mathit{j},\mathit{i}+1},$

$\mathit{M}{{\mathit{ }}^{\mathrm{\prime \prime }}}_{{\mathrm{\Delta }}_j}(t_i)=M_{j,i}^{\ast },\mathit{M}{{\mathit{ }}^{\mathrm{\prime \prime }}}_{{\mathrm{\Delta }}_j}(t_{\mathit{i}+1})=M_{j,i+1}^{\ast },$

to determine $a_{\mathit{j},\mathit{i}},b_{\mathit{j},\mathit{i}},c_{\mathit{j},\mathit{i}}$ and $d_{\mathit{j},\mathit{i}}.$

We get $a_{\mathit{j},\mathit{i}},b_{\mathit{j},\mathit{i}},c_{\mathit{j},\mathit{i}}$ and $d_{\mathit{j},\mathit{i}}$ in (3) as,

(4) $\left.\begin{array}{c}{a}_{\mathit{j},\mathit{i}}=\frac{{\omega }^2{Z}_{\mathit{j},\mathit{i}+1}-M_{j,i+1}^{\ast }+M_{j,i}^{\ast }cos\mathit{\phi }-Z_{\mathit{j},\mathit{i}}{\omega }^{2}cos\mathit{\phi }}{(2sin\mathit{\phi }){\omega }^2},\\ b_{\mathit{j},\mathit{i}}=\frac{{\mathit{Z}}_{\mathit{j},\mathit{i}}{\omega }^2-M_{j,i}^{\ast }}{2{\omega }^2},\\ c_{\mathit{j},\mathit{i}}=\frac{\left({\omega }^2{Z}_{\mathit{j},\mathit{i}+1}+M_{j,i}^{\ast }\right)\left(e^{\phi }+e^{-\mathit{\phi }}\right)-2Z_{\mathit{j},\mathit{i}+1}{\omega }^2-2M_{j,i+1}^{\ast }}{4(e^{-\mathit{\phi }}-e^{\phi }){\omega }^2},\\ d_{\mathit{j},\mathit{i}}=\frac{{\mathit{Z}}_{\mathit{j},\mathit{i}}{\omega }^2+M_{j,i}^{\ast }}{4{\omega }^2},\end{array}\right\}$(4)

where $\mathit{\phi }=\mathit{\omega }(\mathit{t}-t_i)$.

Applying the first derivative continuity condition at $t_i$, $\mathit{M}{{\mathit{ }}^{\mathrm{\prime }}}_{{\mathrm{\Delta }}_j^{-}}(t_i)=\mathit{M}{{\mathit{ }}^{\mathrm{\prime }}}_{{\mathrm{\Delta }}_j^{+}}(t_i)$, we get,

(5) $a_{\mathit{j},\mathit{i}-1}cos\mathit{\phi }-a_{\mathit{j},\mathit{i}}+c_{\mathit{j},\mathit{i}-1}\left(e^{\phi }-e^{-\mathit{\phi }}\right)-2c_{\mathit{j},\mathit{i}}=-b_{\mathit{j},\mathit{i}-1}sin\mathit{\phi }-d_{\mathit{j},\mathit{i}-1}\left(e^{\phi }-e^{-\mathit{\phi }}\right).$(5)

We obtain the following equation by reducing indices in (4) by one and putting them into (5),

(6) $\begin{array}{rl}& \left(\frac{1}{{\mathit{e}}^{-\mathit{\phi }}-e^{\phi }}-\frac{1}{2sin\mathit{\phi }}\right)Z_{\mathit{j},\mathit{i}-1}+\left(cot\mathit{\phi }+\frac{e^{\phi }+e^{-\mathit{\phi }}}{e^{\phi }-e^{-\mathit{\phi }}}\right)Z_{\mathit{j},\mathit{i}}+\left(\frac{1}{{\mathit{e}}^{-\mathit{\phi }}-e^{\phi }}-\frac{1}{2sin\mathit{\phi }}\right)Z_{\mathit{j},\mathit{i}+1}\\ & =\left(\frac{h^2}{{\phi }^2\left(e^{-\mathit{\phi }}-e^{\phi }\right)}-\frac{h^2}{2{\phi }^{2}sin\mathit{\phi }}\right)M_{j,i-1}^{\ast }+\left(\frac{h^{2}cot\mathit{\phi }}{{\phi }^2}+\frac{h^2\left(e^{\phi }+e^{-\mathit{\phi }}\right)}{{\phi }^2\left(e^{\phi }-e^{-\mathit{\phi }}\right)}\right)M_{j,i}^{\ast }\\ & +\left(\frac{h^2}{{\phi }^2\left(e^{-\mathit{\phi }}-e^{\phi }\right)}-\frac{h^2}{2{\phi }^{2}sin\mathit{\phi }}\right)M_{j,i+1}^{\ast }.\end{array}$(6)

Multiplying with $\frac{2\left(e^{\phi }-e^{-\mathit{\phi }}\right)sin\mathit{\phi }}{e^{\phi }-e^{-\mathit{\phi }}+2sin\mathit{\phi }}$ in both sides of (6), we get,

(7) $Z_{\mathit{j},\mathit{i}-1}+\mathit{\rho }{Z}_{\mathit{j},\mathit{i}}+Z_{\mathit{j},\mathit{i}+1}=h^2\left({\alpha }_1{M}_{j,i-1}^{\ast }+{\alpha }_2{M}_{j,i}^{\ast }+{\alpha }_1{M}_{j,i+1}^{\ast }\right),$(7)

where,

$\mathit{\rho }=\frac{2\left(cos\mathit{\phi }+sin\mathit{\phi }\right)e^{\phi }+2\left(sin\mathit{\phi }-cos\mathit{\phi }\right)e^{-\mathit{\phi }}}{e^{-\mathit{\phi }}-e^{\phi }-2sin\mathit{\phi }},$

${\alpha }_1=\frac{e^{\phi }-e^{-\mathit{\phi }}-2sin\mathit{\phi }}{\left(e^{\phi }-e^{-\mathit{\phi }}+2sin\mathit{\phi }\right){\phi }^2},$

${\alpha }_2=\frac{2\left(cos\mathit{\phi }-sin\mathit{\phi }\right)e^{\phi }-2\left(sin\mathit{\phi }+cos\mathit{\phi }\right)e^{-\mathit{\phi }}}{\left(e^{-\mathit{\phi }}-e^{\phi }-2sin\mathit{\phi }\right){\phi }^2}.$

As $\mathit{h}\to 0,\mathit{\phi }\to 0$. We get $\left({\alpha }_1,{\alpha }_2,\mathit{\rho }\right)\to \left(\frac{1}{6},\frac{4}{6},-2\right)$ by applying L’Hospital’s rule. Substituting $\mathit{M}{{\mathit{ }}^{\mathrm{\prime \prime }}}_{{\mathrm{\Delta }}_j}(t_i)=\mathit{Z}{{\mathit{ }}^{\mathrm{\prime \prime }}}_{\mathit{j},\mathit{i}}=M_{j,i}^{\ast }$ into (1). The boundary conditions can be written as,

$Z_{\mathit{j},\mathit{i}}={\varphi }_{\mathit{j},\mathit{i}},-\mathit{N}\le \mathit{i}\le 0,$

$Z_{\mathit{j},2\mathit{N}}={\gamma }_j,$

where $\mathit{j}=1,2$.

For $\mathit{i}=1,2,\dots ,2\mathit{N}-1,\mathit{k}=1,2$, we consider the notations,

$\begin{array}{cccccc}{p}_1(t_i)& =& p_{1,\mathit{i}},& p_2(t_i)& =& p_{2,\mathit{i}},\\ q_{1\mathit{k}}(t_i)& =& q_{1\mathit{k},\mathit{i}},& q_{2\mathit{k}}(t_i)& =& q_{2\mathit{k},\mathit{i}},\\ c_{1\mathit{k}}(t_i)& =& c_{1\mathit{k},\mathit{i}},& c_{2\mathit{k}}(t_i)& =& c_{2\mathit{k},\mathit{i}},\\ s_1(t_i)& =& s_{1,\mathit{i}},& s_2(t_i)& =& s_{2,\mathit{i}}.\end{array}$

From equation (1)(1) $\left.\begin{array}{c}{\mathcal{L}}_1\mathbf{z}\mathbf{(}\mathbf{t}\mathbf{)}\equiv -\mathit{\epsilon }{\mathit{z}{\mathit{ }}^{\mathrm{\prime \prime }}}_1(\mathit{t})+p_1(\mathit{t}){\mathit{z}{\mathit{ }}^{\mathrm{\prime }}}_1(\mathit{t})+{\sum }_{k=1}^2{q}_{1\mathit{k}}(\mathit{t})z_k(\mathit{t})+{\sum }_{k=1}^2{r}_{1\mathit{k}}(\mathit{t})z_k(\mathit{t}-1)=s_1(\mathit{t}),\mathit{t}\in (0,2)\\ {\mathcal{L}}_2\mathbf{z}\mathbf{(}\mathbf{t}\mathbf{)}\equiv -\mathit{\epsilon }{\mathit{z}{\mathit{ }}^{\mathrm{\prime \prime }}}_2(\mathit{t})+p_2(\mathit{t}){\mathit{z}{\mathit{ }}^{\mathrm{\prime }}}_2(\mathit{t})+{\sum }_{k=1}^2{q}_{2\mathit{k}}(\mathit{t})z_k(\mathit{t})+{\sum }_{k=1}^2{r}_{2\mathit{k}}(\mathit{t})z_k(\mathit{t}-1)=s_2(\mathit{t}),\mathit{t}\in (0,2)\end{array}\right\}$(1) , we obtain,

$\mathit{\epsilon }{M}_{1,k}^{\ast }=p_{1,\mathit{k}}\mathit{Z}{{\mathit{ }}^{\mathrm{\prime }}}_{1,\mathit{k}}+q_{11,\mathit{k}}{Z}_{1,\mathit{k}}+q_{12,\mathit{k}}{Z}_{2,\mathit{k}}+r_{11,\mathit{k}}{Z}_1(t_k-1)+r_{12,\mathit{k}}{Z}_2(t_k-1)-s_{1,\mathit{k}},$

$\mathit{\epsilon }{M}_{2,k}^{\ast }=p_{2,\mathit{k}}\mathit{Z}{{\mathit{ }}^{\mathrm{\prime }}}_{2,\mathit{k}}+q_{21,\mathit{k}}{Z}_{2,\mathit{k}}+q_{22,\mathit{k}}{Z}_{2,\mathit{k}}+r_{21,\mathit{k}}{Z}_1(t_k-1)+r_{22,\mathit{k}}{Z}_2(t_k-1)-s_{2,\mathit{k}}.$

Substituting $M_{1,k}^{\ast }$ and $M_{2,k}^{\ast }$ with $\mathit{k}=1,1+\mathit{i},1-\mathit{i}$ and using the following approximations for the first derivative of $\mathit{Z}$,

$\mathit{Z}{{\mathit{ }}^{\mathrm{\prime }}}_{\mathit{j},\mathit{i}}=\frac{{\mathit{Z}}_{\mathit{j},\mathit{i}+1}-Z_{\mathit{j},\mathit{i}-1}}{2\mathit{h}},\mathit{j}=1,2,$

$\mathit{Z}{{\mathit{ }}^{\mathrm{\prime }}}_{\mathit{j},\mathit{i}+1}=\frac{3{\mathit{Z}}_{\mathit{j},\mathit{i}+1}-4Z_{\mathit{j},\mathit{i}}+Z_{\mathit{j},\mathit{i}-1}}{2\mathit{h}},\mathit{j}=1,2,$

$\mathit{Z}{{\mathit{ }}^{\mathrm{\prime }}}_{\mathit{j},\mathit{i}-1}=\frac{-{\mathit{Z}}_{\mathit{j},\mathit{i}+1}+4Z_{\mathit{j},\mathit{i}}-3Z_{\mathit{j},\mathit{i}-1}}{2\mathit{h}},\mathit{j}=1,2.$

We obtain the following linear system of equations for $Z_{1,\mathit{i}}$ and $Z_{2,\mathit{i}},$

$\left.\begin{array}{c}\left[-\mathit{\epsilon }-\frac{3{\alpha }_1\mathit{h}{p}_{1,\mathit{i}-1}}{2}+{\alpha }_1{h}^2{q}_{11,\mathit{i}-1}-\frac{{\alpha }_2\mathit{h}{p}_{1,\mathit{i}}}{2}+\frac{{\alpha }_1\mathit{h}{p}_{1,\mathit{i}+1}}{2}\right]Z_{1,\mathit{i}-1}+\\ \\ \left[-\mathit{\epsilon \rho }+2{\alpha }_1\mathit{h}{p}_{1,\mathit{i}-1}+{\alpha }_2{h}^2{q}_{11,\mathit{i}}-2{\alpha }_1\mathit{h}{p}_{1,\mathit{i}+1}\right]Z_{1,\mathit{i}}\\ \\ +\left[-\mathit{\epsilon }-\frac{{\alpha }_1\mathit{h}{p}_{1,\mathit{i}-1}}{2}+\frac{{\alpha }_2\mathit{h}{p}_{1,\mathit{i}}}{2}+\frac{3{\alpha }_1\mathit{h}{p}_{1,\mathit{i}+1}}{2}+{\alpha }_1{h}^2{q}_{11,\mathit{i}+1}\right]Z_{1,\mathit{i}+1}\\ \\ +{\alpha }_1{h}^2{q}_{12,\mathit{i}-1}{Z}_{2,\mathit{i}-1}+{\alpha }_2{h}^2{q}_{12,\mathit{i}}{Z}_{2,\mathit{i}}+{\alpha }_1{h}^2{q}_{12,\mathit{i}+1}{Z}_{2,\mathit{i}+1}\\ \\ =h^2\left[{\alpha }_1{s}_{1,\mathit{i}-1}+{\alpha }_2{s}_{1,\mathit{i}}+{\alpha }_1{s}_{1,\mathit{i}+1}\right]\\ \\ -h^2\left[{\alpha }_1{r}_{11,\mathit{i}-1}{Z}_1(t_{\mathit{i}-1}-1)+{\alpha }_2{r}_{11,\mathit{i}}{Z}_1(t_i-1)+{\alpha }_1{r}_{11,\mathit{i}+1}{Z}_1(t_{\mathit{i}+1}-1)\right]\\ \\ -h^2\left[{\alpha }_1{r}_{12,\mathit{i}-1}{Z}_2(t_{\mathit{i}-1}-1)+{\alpha }_2{r}_{12,\mathit{i}}{Z}_2(t_i-1)+{\alpha }_1{r}_{12,\mathit{i}+1}{Z}_2(t_{\mathit{i}+1}-1)\right],\end{array}\right\}$

$\left.\begin{array}{c}\left[-\mathit{\epsilon }-\frac{3{\alpha }_1\mathit{h}{p}_{2,\mathit{i}-1}}{2}+{\alpha }_1{h}^2{q}_{22,\mathit{i}-1}-\frac{{\alpha }_2\mathit{h}{p}_{2,\mathit{i}}}{2}+\frac{{\alpha }_1\mathit{h}{p}_{2,\mathit{i}+1}}{2}\right]Z_{2,\mathit{i}-1}+\\ \\ \left[-\mathit{\epsilon \rho }+2{\alpha }_1\mathit{h}{p}_{2,\mathit{i}-1}+{\alpha }_2{h}^2{q}_{22,\mathit{i}}-2{\alpha }_1\mathit{h}{p}_{2,\mathit{i}+1}\right]Z_{2,\mathit{i}}\\ \\ +\left[-\mathit{\epsilon }-\frac{{\alpha }_1\mathit{h}{p}_{2,\mathit{i}-1}}{2}+\frac{{\alpha }_2\mathit{h}{p}_{2,\mathit{i}}}{2}+\frac{3{\alpha }_1\mathit{h}{p}_{2,\mathit{i}+1}}{2}+{\alpha }_1{h}^2{q}_{22,\mathit{i}+1}\right]Z_{2,\mathit{i}+1}\\ \\ +{\alpha }_1{h}^2{q}_{21,\mathit{i}-1}{Z}_{1,\mathit{i}-1}+{\alpha }_2{h}^2{q}_{21,\mathit{i}}{Z}_{1,\mathit{i}}+{\alpha }_1{h}^2{q}_{21,\mathit{i}+1}{Z}_{1,\mathit{i}+1}\\ \\ =h^2\left[{\alpha }_1{s}_{2,\mathit{i}-1}+{\alpha }_2{s}_{2,\mathit{i}}+{\alpha }_1{s}_{2,\mathit{i}+1}\right]\\ \\ -h^2\left[{\alpha }_1{r}_{22,\mathit{i}-1}{Z}_2(t_{\mathit{i}-1}-1)+{\alpha }_2{r}_{22,\mathit{i}}{Z}_2(t_i-1)+{\alpha }_1{r}_{22,\mathit{i}+1}{Z}_2(t_{\mathit{i}+1}-1)\right]\\ \\ -h^2\left[{\alpha }_1{r}_{21,\mathit{i}-1}{Z}_1(t_{\mathit{i}-1}-1)+{\alpha }_2{r}_{21,\mathit{i}}{Z}_1(t_i-1)+{\alpha }_1{r}_{21,\mathit{i}+1}{Z}_1(t_{\mathit{i}+1}-1)\right],\end{array}\right\}$

(8) $\mathrm{for}\mathit{i}=1,2,\dots ,2\mathit{N}-1.$(8)

4.3. Calculation of fitting parameter

An estimate for the solution of the homogeneous problem (1) is of the form:

(9) $z_j(\mathit{t})=z_{j_0}(\mathit{t})+\frac{p_j(0)}{p_j(\mathit{t})}\left(\mathit{\alpha }-z_{j_0}(0)\right)e^{-{\int }_0^t\left(\frac{p_j(\mathit{t})}{\mathit{ϵ}}\right)\mathit{dt}}+\mathit{O}(\mathit{\epsilon }),$(9)

$\mathit{j}=1,2,$ which is derived from the theory of singular perturbations. Here $z_{j_0}(\mathit{t})$ is the solution of the reduced problem.

We may assume that the coefficients are locally constant because we are examining the differential equations on suitably small subintervals. By applying the Taylor’s series expansion for $p_j(\mathit{t}),b_{\mathit{j}1}(\mathit{t})$ and $b_{\mathit{j}2}(\mathit{t})$ about the point $\mathit{t}=0$ and restricting to their first terms, the above equation (9)(9) $z_j(\mathit{t})=z_{j_0}(\mathit{t})+\frac{p_j(0)}{p_j(\mathit{t})}\left(\mathit{\alpha }-z_{j_0}(0)\right)e^{-{\int }_0^t\left(\frac{p_j(\mathit{t})}{\mathit{ϵ}}\right)\mathit{dt}}+\mathit{O}(\mathit{\epsilon }),$(9) becomes,

(10) $z_j(\mathit{t})=z_{j_0}(\mathit{t})+\left(\mathit{\alpha }-z_{j_0}(0)\right)e^{-\left(\frac{p_j(0)}{\mathit{\epsilon }}\right)\mathit{x}}+\mathit{O}(\mathit{\epsilon }),\mathit{j}=1,2.$(10)

From (10), taking the limit as $\mathit{h}\to 0$, and $\mathit{\rho }=\frac{h}{\epsilon }$, we obtain,

$\underset{\mathit{h}\to 0}{lim}{z}_j(\mathit{ih})=z_{j_0}(0)+\left(\mathit{\alpha }-z_{j_0}(0)\right)e^{-\left(p_j(0)\right)\mathit{i\rho }},$

$\underset{\mathit{h}\to 0}{lim}{z}_j(\mathit{ih}+\mathit{h})=z_{j_0}(0)+\left(\mathit{\alpha }-z_{j_0}(0)\right)e^{-\left(p_j(0)\right)(\mathit{i\rho }+\mathit{\rho })},$

$\underset{\mathit{h}\to 0}{lim}{z}_j(\mathit{ih}-\mathit{h})=z_{j_0}(0)+\left(\mathit{\alpha }-z_{j_0}(0)\right)e^{-\left(p_j(0)\right)(\mathit{i\rho }-\mathit{\rho })}.$

Substituting these limiting values in (1), we get the fitting parameter,

${\sigma }_j=\frac{p_j(\mathit{t})\mathit{\rho }}{2}coth\left(\frac{p_j(\mathit{t})\mathit{\rho }}{2}\right),\mathit{j}=1,2.$

The scheme in (8) with the fitting factor can be rewritten as:

$\left.\begin{array}{c}\left[-\mathit{\epsilon }{\sigma }_1-\frac{3{\alpha }_1\mathit{h}{p}_{1,\mathit{i}-1}}{2}+{\alpha }_1{h}^2{q}_{11,\mathit{i}-1}-\frac{{\alpha }_2\mathit{h}{p}_{1,\mathit{i}}}{2}+\frac{{\alpha }_1\mathit{h}{p}_{1,\mathit{i}+1}}{2}\right]Z_{1,\mathit{i}-1}+\\ \\ \left[-\mathit{\epsilon \rho }{\sigma }_1+2{\alpha }_1\mathit{h}{p}_{1,\mathit{i}-1}+{\alpha }_2{h}^2{q}_{11,\mathit{i}}-2{\alpha }_1\mathit{h}{p}_{1,\mathit{i}+1}\right]Z_{1,\mathit{i}}\\ \\ +\left[-\mathit{\epsilon }{\sigma }_1-\frac{{\alpha }_1\mathit{h}{p}_{1,\mathit{i}-1}}{2}+\frac{{\alpha }_2\mathit{h}{p}_{1,\mathit{i}}}{2}+\frac{3{\alpha }_1\mathit{h}{p}_{1,\mathit{i}+1}}{2}+{\alpha }_1{h}^2{q}_{11,\mathit{i}+1}\right]Z_{1,\mathit{i}+1}\\ \\ +{\alpha }_1{h}^2{q}_{12,\mathit{i}-1}{Z}_{2,\mathit{i}-1}+{\alpha }_2{h}^2{q}_{12,\mathit{i}}{Z}_{2,\mathit{i}}+{\alpha }_1{h}^2{q}_{12,\mathit{i}+1}{Z}_{2,\mathit{i}+1}\\ \\ =h^2\left[{\alpha }_1{s}_{1,\mathit{i}-1}+{\alpha }_2{s}_{1,\mathit{i}}+{\alpha }_1{s}_{1,\mathit{i}+1}\right]\\ \\ -h^2\left[{\alpha }_1{r}_{11,\mathit{i}-1}{Z}_1(t_{\mathit{i}-1}-1)+{\alpha }_2{r}_{11,\mathit{i}}{Z}_1(t_i-1)+{\alpha }_1{r}_{11,\mathit{i}+1}{Z}_1(t_{\mathit{i}+1}-1)\right]\\ \\ -h^2\left[{\alpha }_1{r}_{12,\mathit{i}-1}{Z}_2(t_{\mathit{i}-1}-1)+{\alpha }_2{r}_{12,\mathit{i}}{Z}_2(t_i-1)+{\alpha }_1{r}_{12,\mathit{i}+1}{Z}_2(t_{\mathit{i}+1}-1)\right],\end{array}\right\}$

$\left.\begin{array}{c}\left[-\mathit{\epsilon }{\sigma }_2-\frac{3{\alpha }_1\mathit{h}{p}_{2,\mathit{i}-1}}{2}+{\alpha }_1{h}^2{q}_{22,\mathit{i}-1}-\frac{{\alpha }_2\mathit{h}{p}_{2,\mathit{i}}}{2}+\frac{{\alpha }_1\mathit{h}{p}_{2,\mathit{i}+1}}{2}\right]Z_{2,\mathit{i}-1}+\\ \\ \left[-\mathit{\epsilon \rho }{\sigma }_2+2{\alpha }_1\mathit{h}{p}_{2,\mathit{i}-1}+{\alpha }_2{h}^2{q}_{22,\mathit{i}}-2{\alpha }_1\mathit{h}{p}_{2,\mathit{i}+1}\right]Z_{2,\mathit{i}}\\ \\ +\left[-\mathit{\epsilon }{\sigma }_2-\frac{{\alpha }_1\mathit{h}{p}_{2,\mathit{i}-1}}{2}+\frac{{\alpha }_2\mathit{h}{p}_{2,\mathit{i}}}{2}+\frac{3{\alpha }_1\mathit{h}{p}_{2,\mathit{i}+1}}{2}+{\alpha }_1{h}^2{q}_{22,\mathit{i}+1}\right]Z_{2,\mathit{i}+1}\\ \\ +{\alpha }_1{h}^2{q}_{21,\mathit{i}-1}{Z}_{1,\mathit{i}-1}+{\alpha }_2{h}^2{q}_{21,\mathit{i}}{Z}_{1,\mathit{i}}+{\alpha }_1{h}^2{q}_{21,\mathit{i}+1}{Z}_{1,\mathit{i}+1}\\ \\ =h^2\left[{\alpha }_1{s}_{2,\mathit{i}-1}+{\alpha }_2{s}_{2,\mathit{i}}+{\alpha }_1{s}_{2,\mathit{i}+1}\right]\\ \\ -h^2\left[{\alpha }_1{r}_{22,\mathit{i}-1}{Z}_2(t_{\mathit{i}-1}-1)+{\alpha }_2{r}_{22,\mathit{i}}{Z}_2(t_i-1)+{\alpha }_1{r}_{22,\mathit{i}+1}{Z}_2(t_{\mathit{i}+1}-1)\right]\\ \\ -h^2\left[{\alpha }_1{r}_{21,\mathit{i}-1}{Z}_1(t_{\mathit{i}-1}-1)+{\alpha }_2{r}_{21,\mathit{i}}{Z}_1(t_i-1)+{\alpha }_1{r}_{21,\mathit{i}+1}{Z}_1(t_{\mathit{i}+1}-1)\right],\end{array}\right\}$

(11) $\mathrm{for}\mathit{i}=1,2,\dots ,2\mathit{N}-1.$(11)

By incorporating history functions in interval conditions, the above scheme (11) can be rewritten as follows:

$\begin{array}{rl}& E_{1,\mathit{i}}{Z}_{1,\mathit{i}-1}+F_{1,\mathit{i}}{Z}_{1,\mathit{i}}+G_{1,\mathit{i}}{Z}_{1,\mathit{i}+1}+A_{2,\mathit{i}}{Z}_{2,\mathit{i}-1}+B_{2,\mathit{i}}{Z}_{2,\mathit{i}}+C_{2,\mathit{i}}{Z}_{2,\mathit{i}+1}\\ & =D_{1,\mathit{i}}-H_{1,\mathit{i}}{\varphi }_{1,\mathit{i}-1-\mathit{N}}-I_{1,\mathit{i}}{\varphi }_{1,\mathit{i}-\mathit{N}}-J_{1,\mathit{i}}{\varphi }_{1,\mathit{i}+1-\mathit{N}}\\ & -K_{2,\mathit{i}}{\varphi }_{2,\mathit{i}-1-\mathit{N}}-L_{2,\mathit{i}}{\varphi }_{2,\mathit{i}-\mathit{N}}-M_{2,\mathit{i}}{\varphi }_{2,\mathit{i}+1-\mathit{N}},\mathrm{for}1\le \mathit{i}\le \mathit{N}-1,\end{array}$

$\begin{array}{rl}& E_{1,\mathit{i}}{Z}_{1,\mathit{i}-1}+F_{1,\mathit{i}}{Z}_{1,\mathit{i}}+G_{1,\mathit{i}}{Z}_{1,\mathit{i}+1}+A_{2,\mathit{i}}{Z}_{2,\mathit{i}-1}+B_{2,\mathit{i}}{Z}_{2,\mathit{i}}+C_{2,\mathit{i}}{Z}_{2,\mathit{i}+1}\\ & +J_{1,\mathit{i}}{Z}_{1,\mathit{i}+1-\mathit{N}}+M_{2,\mathit{i}}{Z}_{2,\mathit{i}+1-\mathit{N}}=D_{1,\mathit{i}}-H_{1,\mathit{i}}{\varphi }_{1,\mathit{i}-1-\mathit{N}}-I_{1,\mathit{i}}{\varphi }_{1,\mathit{i}-\mathit{N}}\\ & -K_{2,\mathit{i}}{\varphi }_{2,\mathit{i}-1-\mathit{N}}-L_{2,\mathit{i}}{\varphi }_{2,\mathit{i}-\mathit{N}},\mathrm{for}\mathit{i}=\mathit{N},\end{array}$

$\begin{array}{rl}& E_{1,\mathit{i}}{Z}_{1,\mathit{i}-1}+F_{1,\mathit{i}}{Z}_{1,\mathit{i}}+G_{1,\mathit{i}}{Z}_{1,\mathit{i}+1}+A_{2,\mathit{i}}{Z}_{2,\mathit{i}-1}+B_{2,\mathit{i}}{Z}_{2,\mathit{i}}+C_{2,\mathit{i}}{Z}_{2,\mathit{i}+1}\\ & +J_{1,\mathit{i}}{Z}_{1,\mathit{i}+1-\mathit{N}}+M_{2,\mathit{i}}{Z}_{2,\mathit{i}+1-\mathit{N}}+I_{1,\mathit{i}}{Z}_{1,\mathit{i}-\mathit{N}}+L_{2,\mathit{i}}{Z}_{2,\mathit{i}-\mathit{N}}\\ & =D_{1,\mathit{i}}-H_{1,\mathit{i}}{\varphi }_{1,\mathit{i}-1-\mathit{N}}-K_{2,\mathit{i}}{\varphi }_{2,\mathit{i}-1-\mathit{N}}\mathrm{for}\mathit{i}=\mathit{N}+1,\end{array}$

$\begin{array}{rl}& E_{1,\mathit{i}}{Z}_{1,\mathit{i}-1}+F_{1,\mathit{i}}{Z}_{1,\mathit{i}}+G_{1,\mathit{i}}{Z}_{1,\mathit{i}+1}+A_{2,\mathit{i}}{Z}_{2,\mathit{i}-1}+B_{2,\mathit{i}}{Z}_{2,\mathit{i}}+C_{2,\mathit{i}}{Z}_{2,\mathit{i}+1}\\ & +J_{1,\mathit{i}}{Z}_{1,\mathit{i}+1-\mathit{N}}+M_{2,\mathit{i}}{Z}_{2,\mathit{i}+1-\mathit{N}}+I_{1,\mathit{i}}{Z}_{1,\mathit{i}-\mathit{N}}+L_{2,\mathit{i}}{Z}_{2,\mathit{i}-\mathit{N}}\\ & +H_{1,\mathit{i}}{Z}_{1,\mathit{i}-1-\mathit{N}}+K_{2,\mathit{i}}{Z}_{2,\mathit{i}-1-\mathit{N}}=D_{1,\mathit{i}}\mathrm{for}\mathit{N}+2\le \mathit{i}\le 2\mathit{N}-1,\end{array}$

and

$\begin{array}{rl}& E_{2,\mathit{i}}{Z}_{2,\mathit{i}-1}+F_{2,\mathit{i}}{Z}_{2,\mathit{i}}+G_{2,\mathit{i}}{Z}_{2,\mathit{i}+1}+A_{1,\mathit{i}}{Z}_{1,\mathit{i}-1}+B_{1,\mathit{i}}{Z}_{1,\mathit{i}}+C_{1,\mathit{i}}{Z}_{1,\mathit{i}+1}\\ & =D_{2,\mathit{i}}-H_{2,\mathit{i}}{\varphi }_{2,\mathit{i}-1-\mathit{N}}-I_{2,\mathit{i}}{\varphi }_{2,\mathit{i}-\mathit{N}}-J_{2,\mathit{i}}{\varphi }_{2,\mathit{i}+1-\mathit{N}}\\ & -K_{1,\mathit{i}}{\varphi }_{1,\mathit{i}-1-\mathit{N}}-L_{1,\mathit{i}}{\varphi }_{1,\mathit{i}-\mathit{N}}-M_{1,\mathit{i}}{\varphi }_{1,\mathit{i}+1-\mathit{N}}\mathrm{for}1\le \mathit{i}\le \mathit{N}-1,\end{array}$

$\begin{array}{rl}& E_{2,\mathit{i}}{Z}_{2,\mathit{i}-1}+F_{2,\mathit{i}}{Z}_{2,\mathit{i}}+G_{2,\mathit{i}}{Z}_{2,\mathit{i}+1}+A_{1,\mathit{i}}{Z}_{1,\mathit{i}-1}+B_{1,\mathit{i}}{Z}_{1,\mathit{i}}+C_{1,\mathit{i}}{Z}_{1,\mathit{i}+1}\\ & +J_{2,\mathit{i}}{Z}_{2,\mathit{i}+1-\mathit{N}}+M_{1,\mathit{i}}{Z}_{1,\mathit{i}+1-\mathit{N}}=D_{2,\mathit{i}}-K_{2,\mathit{i}}{\varphi }_{2,\mathit{i}-1-\mathit{N}}-I_{2,\mathit{i}}{\varphi }_{2,\mathit{i}-\mathit{N}}\\ & -K_{1,\mathit{i}}{\varphi }_{1,\mathit{i}-1-\mathit{N}}-L_{1,\mathit{i}}{\varphi }_{1,\mathit{i}-\mathit{N}}\mathrm{for}\mathit{i}=\mathit{N},\end{array}$

$\begin{array}{rl}& E_{2,\mathit{i}}{Z}_{2,\mathit{i}-1}+F_{2,\mathit{i}}{Z}_{2,\mathit{i}}+G_{2,\mathit{i}}{Z}_{2,\mathit{i}+1}+A_{1,\mathit{i}}{Z}_{1,\mathit{i}-1}+B_{1,\mathit{i}}{Z}_{1,\mathit{i}}+C_{1,\mathit{i}}{Z}_{1,\mathit{i}+1}\\ & +J_{2,\mathit{i}}{Z}_{2,\mathit{i}+1-\mathit{N}}+M_{1,\mathit{i}}{Z}_{1,\mathit{i}+1-\mathit{N}}+I_{2,\mathit{i}}{Z}_{2,\mathit{i}-\mathit{N}}+L_{1,\mathit{i}}{Z}_{1,\mathit{i}-\mathit{N}}\\ & =D_{2,\mathit{i}}-H_{2,\mathit{i}}{\varphi }_{2,\mathit{i}-1-\mathit{N}}-K_{1,\mathit{i}}{\varphi }_{1,\mathit{i}-1-\mathit{N}}\mathrm{for}\mathit{i}=\mathit{N}+1,\end{array}$

$\begin{array}{rl}& E_{2,\mathit{i}}{Z}_{2,\mathit{i}-1}+F_{2,\mathit{i}}{Z}_{2,\mathit{i}}+G_{2,\mathit{i}}{Z}_{2,\mathit{i}+1}+A_{1,\mathit{i}}{Z}_{1,\mathit{i}-1}+B_{1,\mathit{i}}{Z}_{1,\mathit{i}}+C_{1,\mathit{i}}{Z}_{1,\mathit{i}+1}\\ & +J_{2,\mathit{i}}{Z}_{2,\mathit{i}+1-\mathit{N}}+M_{1,\mathit{i}}{Z}_{1,\mathit{i}+1-\mathit{N}}+I_{2,\mathit{i}}{Z}_{2,\mathit{i}-\mathit{N}}+L_{1,\mathit{i}}{Z}_{1,\mathit{i}-\mathit{N}}\\ & +H_{2,\mathit{i}}{Z}_{2,\mathit{i}-1-\mathit{N}}+K_{1,\mathit{i}}{Z}_{1,\mathit{i}-1-\mathit{N}}=D_{2,\mathit{i}}\mathrm{for}\mathit{N}+2\le \mathit{i}\le 2\mathit{N}-1,\end{array}$

where,

$E_{1,\mathit{i}}=-\mathit{\epsilon }{\sigma }_1-\frac{3{\alpha }_1\mathit{h}{p}_{1,\mathit{i}-1}}{2}+{\alpha }_1{h}^2{q}_{11,\mathit{i}-1}-\frac{{\alpha }_2\mathit{h}{p}_{1,\mathit{i}}}{2}+\frac{{\alpha }_1\mathit{h}{p}_{1,\mathit{i}+1}}{2},$

$F_{1,\mathit{i}}=-\mathit{\epsilon \rho }{\sigma }_1+2{\alpha }_1\mathit{h}{p}_{1,\mathit{i}-1}+{\alpha }_2{h}^2{q}_{11,\mathit{i}}-2{\alpha }_1\mathit{h}{p}_{1,\mathit{i}+1},$

$G_{1,\mathit{i}}=-\mathit{\epsilon }{\sigma }_1-\frac{{\alpha }_1\mathit{h}{p}_{1,\mathit{i}-1}}{2}+\frac{{\alpha }_2\mathit{h}{p}_{1,\mathit{i}}}{2}+\frac{3{\alpha }_1\mathit{h}{p}_{1,\mathit{i}+1}}{2}+{\alpha }_1{h}^2{q}_{11,\mathit{i}+1},$

$D_{1,\mathit{i}}=h^2\left[{\alpha }_1{s}_{1,\mathit{i}-1}+{\alpha }_2{s}_{1,\mathit{i}}+{\alpha }_1{s}_{1,\mathit{i}+1}\right],$

$\begin{array}{ccc}{A}_{1,\mathit{i}}=h^2{\alpha }_1{q}_{12,\mathit{i}-1},& B_{1,\mathit{i}}=h^2{\alpha }_2{q}_{12,\mathit{i}},& C_{1,\mathit{i}}=h^2{\alpha }_1{q}_{12,\mathit{i}+1},\\ \\ H_{1,\mathit{i}}=h^2{\alpha }_1{r}_{11,\mathit{i}-1},& I_{1,\mathit{i}}=h^2{\alpha }_2{r}_{11,\mathit{i}},& J_{1,\mathit{i}}=h^2{\alpha }_1{r}_{11,\mathit{i}+1},\\ \\ K_{1,\mathit{i}}=h^2{\alpha }_1{r}_{12,\mathit{i}-1},& L_{1,\mathit{i}}=h^2{\alpha }_2{r}_{12,\mathit{i}},& M_{1,\mathit{i}}=h^2{\alpha }_1{r}_{12,\mathit{i}+1},\end{array}$

and,

$E_{2,\mathit{i}}=-\mathit{\epsilon }{\sigma }_2-\frac{3{\alpha }_1\mathit{h}{p}_{2,\mathit{i}-1}}{2}+{\alpha }_1{h}^2{q}_{22,\mathit{i}-1}-\frac{{\alpha }_2\mathit{h}{p}_{2,\mathit{i}}}{2}+\frac{{\alpha }_1\mathit{h}{p}_{2,\mathit{i}+1}}{2},$

$F_{2,\mathit{i}}=-\mathit{\epsilon \rho }{\sigma }_2+2{\alpha }_1\mathit{h}{p}_{2,\mathit{i}-1}+{\alpha }_2{h}^2{q}_{22,\mathit{i}}-2{\alpha }_1\mathit{h}{p}_{2,\mathit{i}+1},$

$G_{2,\mathit{i}}=-\mathit{\epsilon }{\sigma }_2-\frac{{\alpha }_1\mathit{h}{p}_{2,\mathit{i}-1}}{2}+\frac{{\alpha }_2\mathit{h}{p}_{2,\mathit{i}}}{2}+\frac{3{\alpha }_1\mathit{h}{p}_{2,\mathit{i}+1}}{2}+{\alpha }_1{h}^2{q}_{22,\mathit{i}+1},$

$D_{2,\mathit{i}}=h^2\left[{\alpha }_1{s}_{2,\mathit{i}-1}+{\alpha }_2{s}_{2,\mathit{i}}+{\alpha }_1{s}_{2,\mathit{i}+1}\right],$

$\begin{array}{ccc}{A}_{2,\mathit{i}}=h^2{\alpha }_1{q}_{21,\mathit{i}-1},& B_{2,\mathit{i}}=h^2{\alpha }_2{q}_{21,\mathit{i}},& C_{2,\mathit{i}}=h^2{\alpha }_1{q}_{21,\mathit{i}+1},\\ \\ H_{2,\mathit{i}}=h^2{\alpha }_1{r}_{22,\mathit{i}-1},& I_{2,\mathit{i}}=h^2{\alpha }_2{r}_{22,\mathit{i}},& J_{2,\mathit{i}}=h^2{\alpha }_1{r}_{22,\mathit{i}+1},\\ \\ K_{2,\mathit{i}}=h^2{\alpha }_1{r}_{21,\mathit{i}-1},& L_{2,\mathit{i}}=h^2{\alpha }_2{r}_{21,\mathit{i}},& M_{2,\mathit{i}}=h^2{\alpha }_1{r}_{21,\mathit{i}+1}.\end{array}$

5. Error and convergence analysis

By Taylor series expansion,

$z_{\mathit{j},\mathit{i}+1}=z_{\mathit{j},\mathit{i}}+\mathit{h}\mathit{z}{{\mathit{ }}^{\mathrm{\prime }}}_{\mathit{j},\mathit{i}}+\frac{h^2}{2}\mathit{z}{{\mathit{ }}^{\mathrm{\prime \prime }}}_{\mathit{j},\mathit{i}}+\frac{h^3}{6}{z}_{j,i}^3+\mathit{O}(h^4),$

$z_{\mathit{j},\mathit{i}-1}=z_{\mathit{j},\mathit{i}}-\mathit{h}\mathit{z}{{\mathit{ }}^{\mathrm{\prime }}}_{\mathit{j},\mathit{i}}+\frac{h^2}{2}\mathit{z}{{\mathit{ }}^{\mathrm{\prime \prime }}}_{\mathit{j},\mathit{i}}-\frac{h^3}{6}{z}_{j,i}^3+\mathit{O}(h^4).$

For $\mathit{j}=1,2$, the truncation error obtained is

${\tau }_{\mathit{j},\mathit{i}}=\frac{h^3}{3}{\alpha }_1{p}_{\mathit{j},\mathit{i}-1}{z}_{j,i}^3-\frac{h^3}{6}{\alpha }_2{p}_{\mathit{j},\mathit{i}}{z}_{j,i}^3+\frac{h^3}{3}{\alpha }_1{p}_{\mathit{j},\mathit{i}+1}{z}_{j,i}^3+\mathit{O}(h^4).$

We combine the equations derived from (11) for the system of two equations to arrive at a matrix form as follows:

(12) $\mathit{XZ}=\mathit{V},$(12)

where $\mathit{X}=(x_{\mathit{ij}})$ is a matrix of order $4\mathit{N}-2$ such that,

$x_{\mathit{i},\mathit{i}}=\left\{\begin{array}{ccc}{F}_{1,\mathit{i}}& ,& \mathit{i}=1\mathrm{to}2\mathit{N}-1\\ F_{2,\mathit{i}}& ,& \mathit{i}=2\mathit{N}\mathrm{to}4\mathit{N}-1;\end{array}\right.$

$x_{\mathit{i},\mathit{i}+1}=\left\{\begin{array}{ccc}{G}_{1,\mathit{i}}& ,& \mathit{i}=1\mathrm{to}2\mathit{N}-1\\ G_{2,\mathit{i}}& ,& \mathit{i}=2\mathit{N}+1\mathrm{to}4\mathit{N}-2;\end{array}\right.$

$x_{\mathit{i},\mathit{i}-1}=\left\{\begin{array}{ccc}{E}_{1,\mathit{i}}& ,& \mathit{i}=2\mathrm{to}2\mathit{N}-1\\ E_{2,\mathit{i}}& ,& \mathit{i}=2\mathit{N}\mathrm{to}4\mathit{N}-2;\end{array}\right.$

$\begin{array}{cc}{x}_{\mathit{i},\mathit{i}+2\mathit{N}}=C_{2,\mathit{i}},\mathit{i}=1\mathrm{to}2\mathit{N}-2;& x_{\mathit{i},\mathit{i}+2\mathit{N}-1}=B_{2,\mathit{i}},\mathit{i}=1\mathrm{to}2\mathit{N}-1;\\ x_{\mathit{i},\mathit{i}+2\mathit{N}-2}=A_{2,\mathit{i}},\mathit{i}=2\mathrm{to}2\mathit{N}-1;& x_{\mathit{i},\mathit{i}-\mathit{N}}=I_{1,\mathit{i}},\mathit{i}=\mathit{N}+1\mathrm{to}2\mathit{N}-1;\\ x_{\mathit{i},\mathit{i}-\mathit{N}+1}=J_{1,\mathit{i}},\mathit{i}=\mathit{N}+1\mathrm{to}2\mathit{N}-1;& x_{\mathit{i},\mathit{i}-\mathit{N}-1}=H_{1,\mathit{i}},\mathit{i}=\mathit{N}+2\mathrm{to}2\mathit{N}-1;\\ x_{\mathit{i},\mathit{i}+\mathit{N}}=M_{2,\mathit{i}},\mathit{i}=\mathit{N}+1\mathrm{to}2\mathit{N}-1;& x_{\mathit{i},\mathit{i}+\mathit{N}-1}=L_{2,\mathit{i}},\mathit{i}=\mathit{N}+1\mathrm{to}2\mathit{N}-1;\\ x_{\mathit{i},\mathit{i}+\mathit{N}-2}=K_{2,\mathit{i}},\mathit{i}=\mathit{N}+2\mathrm{to}2\mathit{N}-1;& x_{\mathit{i},\mathit{i}-2\mathit{N}}=A_{1,\mathit{i}},\mathit{i}=2\mathit{N}+1\mathrm{to}4\mathit{N}-2;\\ x_{\mathit{i},\mathit{i}-2\mathit{N}+1}=B_{1,\mathit{i}},\mathit{i}=2\mathit{N}\mathrm{to}4\mathit{N}-2;& x_{\mathit{i},\mathit{i}-2\mathit{N}+2}=C_{1,\mathit{i}},\mathit{i}=2\mathit{N}\mathrm{to}4\mathit{N}-2;\\ x_{\mathit{i},\mathit{i}-2\mathit{N}-2}=M_{1,\mathit{i}},\mathit{i}=2\mathit{N}+3\mathrm{to}4\mathit{N}-2;& x_{\mathit{i},\mathit{i}-2\mathit{N}-3}=L_{1,\mathit{i}},\mathit{i}=2\mathit{N}+4\mathrm{to}4\mathit{N}-2;\\ x_{\mathit{i},\mathit{i}-2\mathit{N}-4}=K_{1,\mathit{i}},\mathit{i}=2\mathit{N}+5\mathrm{to}4\mathit{N}-2;& x_{\mathit{i},\mathit{i}-3}=J_{2,\mathit{i}},\mathit{i}=2\mathit{N}+3\mathrm{to}4\mathit{N}-2;\\ x_{\mathit{i},\mathit{i}-4}=I_{2,\mathit{i}},\mathit{i}=2\mathit{N}+4\mathrm{to}4\mathit{N}-2;& x_{\mathit{i},\mathit{i}-5}=H_{2,\mathit{i}},\mathit{i}=2\mathit{N}+5\mathrm{to}4\mathit{N}-2,\end{array}$

and all the remaining entries of $\mathit{X}$ are zero and $\mathit{V}=(v_i)$ is a column vector such that

$v_i=\left\{\begin{array}{cc}{D}_{1,\mathit{i}}-H_{1,\mathit{i}}{\varphi }_{1,\mathit{i}-1-\mathit{N}}-I_{1,\mathit{i}}{\varphi }_{1,\mathit{i}-\mathit{N}}-J_{1,\mathit{i}}{\varphi }_{1,\mathit{i}+1-\mathit{N}}& \mathit{\hspace{0.17em}}\\ -K_{2,\mathit{i}}{\varphi }_{2,\mathit{i}-1-\mathit{N}}-L_{2,\mathit{i}}{\varphi }_{2,\mathit{i}-\mathit{N}}-M_{2,\mathit{i}}{\varphi }_{2,\mathit{i}+1-\mathit{N}},& \mathrm{for}1\le \mathit{i}\le \mathit{N}-1,\\ D_{1,\mathit{i}}-H_{1,\mathit{i}}{\varphi }_{1,\mathit{i}-1-\mathit{N}}-I_{1,\mathit{i}}{\varphi }_{1,\mathit{i}-\mathit{N}}-K_{2,\mathit{i}}{\varphi }_{2,\mathit{i}-1-\mathit{N}}& \mathit{\hspace{0.17em}}\\ -L_{2,\mathit{i}}{\varphi }_{2,\mathit{i}-\mathit{N}},& \mathrm{for}\mathit{i}=\mathit{N},\\ D_{1,\mathit{i}}-H_{1,\mathit{i}}{\varphi }_{1,\mathit{i}-1-\mathit{N}}-K_{2,\mathit{i}}{\varphi }_{2,\mathit{i}-1-\mathit{N}}& \mathrm{for}\mathit{i}=\mathit{N}+1,\\ D_{1,\mathit{i}}& \mathrm{for}\mathit{N}+2\le \mathit{i}\le 2\mathit{N}-1,\\ D_{2,\mathit{i}}-H_{2,\mathit{i}}{\varphi }_{2,\mathit{i}-1-\mathit{N}}-I_{2,\mathit{i}}{\varphi }_{2,\mathit{i}-\mathit{N}}-J_{2,\mathit{i}}{\varphi }_{2,\mathit{i}+1-\mathit{N}}& \mathit{\hspace{0.17em}}\\ -K_{1,\mathit{i}}{\varphi }_{1,\mathit{i}-1-\mathit{N}}-L_{1,\mathit{i}}{\varphi }_{1,\mathit{i}-\mathit{N}}-M_{1,\mathit{i}}{\varphi }_{1,\mathit{i}+1-\mathit{N}}& \mathrm{for}2\mathit{N}\le \mathit{i}\le 3\mathit{N}-1,\\ D_{2,\mathit{i}}-K_{2,\mathit{i}}{\varphi }_{2,\mathit{i}-1-\mathit{N}}-I_{2,\mathit{i}}{\varphi }_{2,\mathit{i}-\mathit{N}}-K_{1,\mathit{i}}{\varphi }_{1,\mathit{i}-1-\mathit{N}}& \mathit{\hspace{0.17em}}\\ -L_{1,\mathit{i}}{\varphi }_{1,\mathit{i}-\mathit{N}},& \mathrm{for}\mathit{i}=3\mathit{N},\\ D_{2,\mathit{i}}-H_{2,\mathit{i}}{\varphi }_{2,\mathit{i}-1-\mathit{N}}-K_{1,\mathit{i}}{\varphi }_{1,\mathit{i}-1-\mathit{N}}& \mathrm{for}\mathit{i}=3\mathit{N}+1,\\ D_{2,\mathit{i}},& \mathrm{for}3\mathit{N}+2\le \mathit{i}\le 4\mathit{N}-2,\end{array}\right.$

$\mathrm{and}\mathit{Z}={\left[\begin{array}{cccccccccc}{\mathit{Z}}_{1,1}(\mathit{h})& .& .& .& Z_{1,2\mathit{N}-1}(\mathit{h})& Z_{2,1}(\mathit{h})& .& .& .& Z_{2,2\mathit{N}-1}(\mathit{h})\end{array}\right]}^{\mathit{t}}.$

Let

$\mathit{Y}(\mathit{h})={\left[\begin{array}{cccccccccc}{\mathit{Y}}_{1,1}(\mathit{h})& .& .& .& Y_{1,2\mathit{N}-1}(\mathit{h})& Y_{2,1}(\mathit{h})& .& .& .& Y_{2,2\mathit{N}-1}(\mathit{h})\end{array}\right]}^{\mathit{t}}$

be the truncation error, where

$Y_{\mathit{j},\mathit{i}}(\mathit{h})=\frac{h^3}{3}\left[{\alpha }_1{p}_{\mathit{j},\mathit{i}-1}{z}_{j,i}^3-{\alpha }_2{p}_{\mathit{j},\mathit{i}}{z}_{j,i}^3+{\alpha }_1{p}_{\mathit{j},\mathit{i}+1}{z}_{j,i}^3\right]+\mathit{O}(h^4).$

Then (12) can be written as $\mathit{X}\bar{Z}-\mathit{Y}(\mathit{h})=\mathit{V}$, where

$\bar{Z}={\left[\begin{array}{cccccccccc}{\bar{Z}}_{1,1}(\mathit{h})& .& .& .& {\bar{Z}}_{1,2\mathit{N}-1}(\mathit{h})& {\bar{Z}}_{2,1}(\mathit{h})& .& .& .& {\bar{Z}}_{2,2\mathit{N}-1}(\mathit{h})\end{array}\right]}^{\mathit{t}}$

is the exact solution. It follows that $\mathit{X}{\bar{E}}_{\mathit{r}}=\mathit{Y}(\mathit{h})$, where

${\bar{E}}_{\mathit{r}}=\bar{Z}-\mathit{Z}={\bar{E}}_{\mathit{r}}={\left[\begin{array}{cccccccccc}{\mathit{e}}_{1,1}& .& .& .& e_{1,2\mathit{N}-1}& e_{2,1}& .& .& .& e_{2,2\mathit{N}-1}\end{array}\right]}^{\mathit{t}}$

Let ${\bar{R}}_{\mathit{i}}$ denote the $i^{\mathit{th}}$ row sum of the matrix $\mathit{X}$. Then

$\begin{array}{ccc}{\bar{R}}_1& =& F_{1,1}+G_{1,1}+B_{2,1}+C_{2,1}\\ \mathit{\hspace{0.17em}}& =& \mathit{\epsilon }{\sigma }_1+\frac{3}{2}{\alpha }_1\mathit{h}{p}_{1,\mathit{i}-1}+\frac{{\alpha }_2\mathit{h}{p}_{1,\mathit{i}}}{2}-\frac{{\alpha }_1\mathit{h}{p}_{1,\mathit{i}+1}}{2}+{\alpha }_2{h}^2{q}_{11,\mathit{i}}+{\alpha }_1{h}^2{q}_{11,\mathit{i}+1}+h^2{\alpha }_2{q}_{21,\mathit{i}}+h^2{\alpha }_1{q}_{21,\mathit{i}+1},\\ \\ {\bar{R}}_{\mathit{i}}& =& E_{1,\mathit{i}}{F}_{1,\mathit{i}}+G_{1,\mathit{i}}+A_{2,\mathit{i}}+B_{2,\mathit{i}}+C_{2,\mathit{i}}\\ \mathit{\hspace{0.17em}}& =& h^2[{\alpha }_1{q}_{11,\mathit{i}-1}+{\alpha }_2{q}_{11,\mathit{i}}+{\alpha }_1{q}_{11,\mathit{i}+1}+{\alpha }_1{q}_{21,\mathit{i}-1}+{\alpha }_2{q}_{21,\mathit{i}}+{\alpha }_1{q}_{21,\mathit{i}+1}]\\ \mathit{\hspace{0.17em}}& =& h^2{Q}_i,\mathrm{for}\mathit{i}=2\mathrm{to}\mathit{N},\\ \\ {\bar{R}}_{\mathit{i}}& =& I_{1,\mathit{i}}+J_{1,\mathit{i}}+E_{1,\mathit{i}}+F_{1,\mathit{i}}+G_{1,\mathit{i}}+L_{2,\mathit{i}}+M_{2,\mathit{i}}+A_{2,\mathit{i}}+B_{2,\mathit{i}}+C_{2,\mathit{i}}\\ \mathit{\hspace{0.17em}}& =& h^2\left[{\alpha }_1{q}_{11,\mathit{i}-1}+{\alpha }_2{q}_{11,\mathit{i}}+{\alpha }_1{q}_{11,\mathit{i}+1}+{\alpha }_1{q}_{21,\mathit{i}-1}+{\alpha }_2{q}_{21,\mathit{i}}\right.\\ \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& \left.+{\alpha }_1{q}_{21,\mathit{i}+1}+{\alpha }_2{r}_{11,\mathit{i}}+{\alpha }_1{r}_{11,\mathit{i}+1}+{\alpha }_2{r}_{21,\mathit{i}}+{\alpha }_1{r}_{21,\mathit{i}+1}\right]\\ \mathit{\hspace{0.17em}}& =& h^2{Q}_i,\mathrm{for}\mathit{i}=\mathit{N}+1,\\ \\ {\bar{R}}_{\mathit{i}}& =& I_{1,\mathit{i}}+J_{1,\mathit{i}}+E_{1,\mathit{i}}+F_{1,\mathit{i}}+G_{1,\mathit{i}}+L_{2,\mathit{i}}+M_{2,\mathit{i}}+A_{2,\mathit{i}}+B_{2,\mathit{i}}+C_{2,\mathit{i}}\\ \mathit{\hspace{0.17em}}& =& h^2\left[{\alpha }_1{q}_{11,\mathit{i}-1}+{\alpha }_2{q}_{11,\mathit{i}}+{\alpha }_1{q}_{11,\mathit{i}+1}+{\alpha }_1{q}_{21,\mathit{i}-1}+{\alpha }_2{q}_{21,\mathit{i}}\right.\\ \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& \left.+{\alpha }_1{q}_{21,\mathit{i}+1}+{\alpha }_2{r}_{11,\mathit{i}}+{\alpha }_1{r}_{11,\mathit{i}+1}+{\alpha }_2{r}_{21,\mathit{i}}+{\alpha }_1{r}_{21,\mathit{i}+1}+{\alpha }_1{r}_{11,\mathit{i}-1}+{\alpha }_2{r}_{21,\mathit{i}-1}\right]\\ \mathit{\hspace{0.17em}}& =& h^2{Q}_i,\mathrm{for}\mathit{i}=\mathit{N}+2\mathrm{to}2\mathit{N}-2,\\ \\ {\bar{R}}_{\mathit{i}}& =& H_{1,\mathit{i}}+I_{1,\mathit{i}}+J_{1,\mathit{i}}+E_{1,\mathit{i}}+F_{1,\mathit{i}}+K_{2,\mathit{i}}+L_{2,\mathit{i}}+M_{2,\mathit{i}}+A_{2,\mathit{i}}+B_{2,\mathit{i}}\\ \mathit{\hspace{0.17em}}& =& \mathit{\epsilon }{\sigma }_1+\frac{{\alpha }_1\mathit{h}{p}_{1,\mathit{i}-1}}{2}+{\alpha }_2{h}^2{q}_{11,\mathit{i}}+{\alpha }_1{h}^2{q}_{11,\mathit{i}-1}-\frac{3}{2}{\alpha }_1\mathit{h}{p}_{1,\mathit{i}+1}-\frac{{\alpha }_2\mathit{h}{p}_{1,\mathit{i}}}{2}+h^2\left[{\alpha }_1{r}_{11,\mathit{i}-1}+{\alpha }_2{r}_{11,\mathit{i}}\right.\\ \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& \left.+{\alpha }_1{r}_{11,\mathit{i}+1}+{\alpha }_1{r}_{21,\mathit{i}-1}+{\alpha }_2{r}_{21,\mathit{i}}+{\alpha }_1{r}_{21,\mathit{i}+1}+{\alpha }_1{q}_{21,\mathit{i}-1}+{\alpha }_2{q}_{21,\mathit{i}}\right],\mathrm{for}\mathit{i}=2\mathit{N}-1,\\ \\ {\bar{R}}_{\mathit{i}}& =& B_{1,\mathit{i}}+C_{1,\mathit{i}}+F_{2,\mathit{i}}+G_{2,\mathit{i}}\\ \mathit{\hspace{0.17em}}& =& h^2[{\alpha }_2{q}_{12,\mathit{i}}+{\alpha }_1{q}_{12,\mathit{i}+1}+{\alpha }_2{q}_{22,\mathit{i}}+{\alpha }_1{q}_{22,\mathit{i}+1}]\\ \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& +\mathit{\epsilon }{\sigma }_2+\frac{3\mathit{h}}{2}{\alpha }_1{p}_{2,\mathit{i}-1}+\frac{h}{2}{\alpha }_2{p}_{2,\mathit{i}}-\frac{h}{2}{\alpha }_1{p}_{2,\mathit{i}+1},\mathrm{for}\mathit{i}=2\mathit{N},\end{array}$

$\begin{array}{ccc}{\bar{R}}_{\mathit{i}}& =& E_{2,\mathit{i}}{F}_{2,\mathit{i}}+G_{2,\mathit{i}}+A_{1,\mathit{i}}+B_{1,\mathit{i}}+C_{1,\mathit{i}}\\ \mathit{\hspace{0.17em}}& =& h^2[{\alpha }_1{q}_{22,\mathit{i}-1}+{\alpha }_2{q}_{22,\mathit{i}}+{\alpha }_1{q}_{22,\mathit{i}+1}+{\alpha }_1{q}_{12,\mathit{i}-1}+{\alpha }_2{q}_{12,\mathit{i}}+{\alpha }_1{q}_{12,\mathit{i}+1}]\\ \mathit{\hspace{0.17em}}& =& h^2{Q}_i,\mathrm{for}\mathit{i}=2\mathit{N}+1\mathrm{to}3\mathit{N}-2,\\ \\ {\bar{R}}_{\mathit{i}}& =& J_{2,\mathit{i}}+E_{2,\mathit{i}}+F_{2,\mathit{i}}+G_{2,\mathit{i}}+M_{1,\mathit{i}}+A_{1,\mathit{i}}+B_{1,\mathit{i}}+C_{1,\mathit{i}}\\ \mathit{\hspace{0.17em}}& =& h^2[{\alpha }_1{q}_{22,\mathit{i}-1}+{\alpha }_2{q}_{22,\mathit{i}}+{\alpha }_1{q}_{22,\mathit{i}+1}+{\alpha }_1{q}_{12,\mathit{i}-1}+{\alpha }_2{q}_{12,\mathit{i}}+{\alpha }_1{q}_{12,\mathit{i}+1}+{\alpha }_1{r}_{12,\mathit{i}+1}+{\alpha }_1{r}_{22,\mathit{i}+1}]\\ \mathit{\hspace{0.17em}}& =& h^2{Q}_i,\mathrm{for}\mathit{i}=3\mathit{N}-1,\\ \\ {\bar{R}}_{\mathit{i}}& =& I_{2,\mathit{i}}+J_{2,\mathit{i}}+E_{2,\mathit{i}}+F_{2,\mathit{i}}+G_{2,\mathit{i}}+L_{1,\mathit{i}}+M_{1,\mathit{i}}+A_{1,\mathit{i}}+B_{1,\mathit{i}}+C_{1,\mathit{i}}\\ \mathit{\hspace{0.17em}}& =& h^2\left[{\alpha }_1{q}_{22,\mathit{i}-1}+{\alpha }_2{q}_{22,\mathit{i}}+{\alpha }_1{q}_{22,\mathit{i}+1}+{\alpha }_1{q}_{12,\mathit{i}-1}+{\alpha }_2{q}_{12,\mathit{i}}+{\alpha }_1{q}_{12,\mathit{i}+1}+{\alpha }_1{r}_{12,\mathit{i}+1}\right.\\ \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& \left.+{\alpha }_1{r}_{22,\mathit{i}+1}+{\alpha }_2{r}_{12,\mathit{i}}+{\alpha }_2{r}_{22,\mathit{i}}\right]\\ \mathit{\hspace{0.17em}}& =& h^2{Q}_i,\mathrm{for}\mathit{i}=3\mathit{N},\\ \\ {\bar{R}}_{\mathit{i}}& =& H_{2,\mathit{i}}+I_{2,\mathit{i}}+J_{2,\mathit{i}}+E_{2,\mathit{i}}+F_{2,\mathit{i}}+G_{2,\mathit{i}}+K_{1,\mathit{i}}+L_{1,\mathit{i}}+M_{1,\mathit{i}}+A_{1,\mathit{i}}+B_{1,\mathit{i}}+C_{1,\mathit{i}}\\ \mathit{\hspace{0.17em}}& =& h^2\left[{\alpha }_1{q}_{22,\mathit{i}-1}+{\alpha }_2{q}_{22,\mathit{i}}+{\alpha }_1{q}_{22,\mathit{i}+1}+{\alpha }_1{q}_{12,\mathit{i}-1}+{\alpha }_2{q}_{12,\mathit{i}}+{\alpha }_1{q}_{12,\mathit{i}+1}+{\alpha }_1{r}_{12,\mathit{i}+1}\right.\\ \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& \left.+{\alpha }_1{r}_{22,\mathit{i}+1}+{\alpha }_2{r}_{12,\mathit{i}}+{\alpha }_2{r}_{[}22,\mathit{i}]+{\alpha }_1{r}_{12,\mathit{i}-1}+{\alpha }_1{r}_{22,\mathit{i}-1}\right]\\ \mathit{\hspace{0.17em}}& =& h^2{Q}_i,\mathrm{for}\mathit{i}=3\mathit{N}+1\mathrm{to}4\mathit{N}-1,\\ \\ {\bar{R}}_{\mathit{i}}& =& H_{2,\mathit{i}}+I_{2,\mathit{i}}+J_{2,\mathit{i}}+E_{2,\mathit{i}}+F_{2,\mathit{i}}+K_{1,\mathit{i}}+L_{1,\mathit{i}}+M_{1,\mathit{i}}+A_{1,\mathit{i}}+B_{1,\mathit{i}}\\ \mathit{\hspace{0.17em}}& =& \mathit{\epsilon }{\sigma }_2+\frac{3{\alpha }_1\mathit{h}{p}_{2,\mathit{i}-1}}{2}-\frac{h}{2}{\alpha }_2{p}_{2,\mathit{i}}-\frac{3\mathit{h}}{2}{\alpha }_1{p}_{2,\mathit{i}}+h^2\left[{\alpha }_1{r}_{12,\mathit{i}-1}+{\alpha }_2{r}_{12,\mathit{i}}+{\alpha }_1{r}_{12,\mathit{i}+1}\right.\\ \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}& \left.+{\alpha }_1{r}_{22,\mathit{i}-1}+{\alpha }_2{r}_{22,\mathit{i}}+{\alpha }_1{r}_{22,\mathit{i}+1}+{\alpha }_1{q}_{22,\mathit{i}-1}+{\alpha }_2{q}_{22,\mathit{i}}+{\alpha }_1{q}_{12,\mathit{i}-1}+{\alpha }_2{q}_{12,\mathit{i}}\right],\mathrm{for}\mathit{i}=4\mathit{N}-2.\end{array}$

It can be easily verified that the matrix $\mathit{X}$ is monotone and irreducible when $\mathit{h}$ is sufficiently small and hence $X^{-1}$ exists and its elements are non-negative. From ${\bar{E}}_{\mathit{r}}=X^{-1}\mathit{Y}(\mathit{h})$ which implies

$||{\bar{E}}_{\mathit{r}}||=||X^{-1}||||\mathit{Y}(\mathit{h})||$

Let ${\bar{x}}_{\mathit{j},\mathit{i}}$ be the $(\mathit{j},\mathit{i}{)}^{\mathit{th}}$ element of $X^{-1}$. Then ${\sum }_{i=1}^{4N-2}{\bar{x}}_{\mathit{j},\mathit{i}}{\bar{R}}_{\mathit{i}}=1$ for $\mathit{j}=1(1)(4\mathit{N}-2)$. This gives

$\sum _{\mathit{i}=1}^{4\mathit{N}-2}{\bar{x}}_{\mathit{j},\mathit{i}}\le \frac{1}{{min}_{1\le \mathit{i}\le 4\mathit{N}-2}{\bar{R}}_{\mathit{i}}}\le \frac{1}{h^2|Q_i|}$

Then we can conclude that, for $\mathit{j}=1,2$ and $\mathit{i}=1,2,\dots ,2\mathit{N}-1$

$e_{\mathit{j},\mathit{i}}\le \frac{1}{h^2|Q_i|}\frac{h^3{L}^{\ast }}{3}=\mathit{O}(\mathit{h})$

Hence $||\bar{E}||=\mathit{O}(\mathit{h})$

6. Numerical experiments and discussions

Some numerical experiments are taken into consideration to illustrate the applicability of the proposed method. The tabulated solution of some problems with varying $\mathit{\epsilon }$ is presented. Our theoretical approaches have been validated by solving the following examples numerically. MATLAB R2022a is utilized to plot the estimated findings. For the problems with no exact solution, maximum absolute errors are determined by applying the double mesh principle to the provided examples.

(13) $E_{i,\epsilon }^N=\underset{0\le \mathit{j}\le \mathit{N}}{max}|Z_{i,j}^N-Z_{i,2j}^{2N}|,$(13)

for $\mathit{i}=1,2,$ where $Z_{i,j}^N$ and $Z_{i,2j}^{2N}$ are the $j^{\mathit{th}}$ and $2j^{\mathit{th}}$ components of the numerical solutions of meshes of $\mathit{N}$ and $2\mathit{N}$ points, respectively. We compute the uniform error by,

$E_i^N=\underset{\epsilon }{max}{E}_{i,\epsilon }^N.$

Example 1

$\left\{\begin{array}{ccc}-\mathit{\epsilon }{\mathit{z}{\mathit{ }}^{\mathrm{\prime \prime }}}_1(\mathit{t})+11z_{1{\mathit{ }}^{\mathrm{\prime }}}(\mathit{t})+6z_1(\mathit{t})-2z_2(\mathit{t})-z_1(\mathit{t}-1)& =& e^t,\mathit{t}\in [0,2],\\ -\mathit{\epsilon }{\mathit{z}{\mathit{ }}^{\mathrm{\prime \prime }}}_2(\mathit{t})+16z_{2{\mathit{ }}^{\mathrm{\prime }}}(\mathit{t})-2z_1(\mathit{t})+5z_2(\mathit{t})-z_2(\mathit{t}-1)& =& t^2,\mathit{t}\in [0,2],\\ z_1(\mathit{t})=1,\mathit{t}\in [-1,0],z_1(2)=1,& \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}\\ z_2(\mathit{t})=1,\mathit{t}\in [-1,0],z_2(2)=1.& \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}\end{array}\right.$

The maximum pointwise errors are presented in Table 1 for different values of $\mathit{\epsilon }$ and the numerical results compared with the results in (Subburayan and Ramanujam 2014; Selvi and Ramanujam 2017; Chakravarthy and Gupta 2020). Table 3 shows the maximum absolute errors for different values of $\mathit{\epsilon }$. As the step size $\mathit{h}=1/\mathit{N}$ tends to zero for all values of the perturbation parameter epsilon, the maximum absolute errors decrease with the linear rate of convergence, as shown by the results presented in Table 2 for different values of $\mathit{N}$ with $\mathit{\epsilon }=0.1$.

Table 3. The maximum absolute error of the example 1 for different values of $\mathit{\epsilon }$.

$\mathit{\epsilon }$	N
	$64$	$128$	$256$	$512$	$1024$	$2048$
$2^{-2}$	2.3734e-03	1.1490e-03	5.5708e-04	2.7285e-04	1.3480e-04	6.6970e-05
$2^{-4}$	2.5442e-03	1.3587e-03	6.7680e-04	3.2860e-04	1.5980e-04	7.8422e-05
$2^{-6}$	1.8228e-03	1.1319e-03	6.6790e-04	3.5628e-04	1.7755e-04	8.6314e-05
$2^{-8}$	1.6444e-03	8.3625e-04	4.6226e-04	2.8751e-04	1.6970e-04	9.0491e-05
$2^{-10}$	1.6443e-03	8.2997e-04	4.1701e-04	2.1059e-04	1.1614e-04	7.2310e-05
$2^{-12}$	1.6443e-03	8.2997e-04	4.1700e-04	2.0901e-04	1.0464e-04	5.2745e-05
$2^{-14}$	1.6443e-03	8.2997e-04	4.1700e-04	2.0901e-04	1.0463e-04	5.2349e-05
$E_1^N$	1.6443e-03	8.2997e-04	4.1700e-04	2.0901e-04	1.0463e-04	5.2349e-05
$2^{-2}$	1.3603e-03	6.2837e-04	2.8960e-04	1.3656e-04	6.5900e-05	3.2311e-05
$2^{-4}$	1.3443e-03	7.4396e-04	3.6648e-04	1.7022e-04	7.8822e-05	3.7283e-05
$2^{-6}$	9.4431e-04	5.6304e-04	3.4369e-04	1.9069e-04	9.4098e-05	4.3789e-05
$2^{-8}$	9.0261e-04	4.5233e-04	2.3649e-04	1.4101e-04	8.6523e-05	4.8060e-05
$2^{-10}$	9.0260e-04	4.5166e-04	2.2593e-04	1.1316e-04	5.9150e-05	3.5268e-05
$2^{-12}$	9.0260e-04	4.5166e-04	2.2593e-04	1.1299e-04	5.6502e-05	2.8295e-05
$2^{-14}$	9.0260e-04	4.5166e-04	2.2593e-04	1.1299e-04	5.6502e-05	2.8253e-05
$E_2^N$	9.0260e-04	4.5166e-04	2.2593e-04	1.1299e-04	5.6502e-05	2.8253e-05

Table 1. The maximum absolute error of the example 1 and results in (Subburayan and Ramanujam 2014; Selvi and Ramanujam 2017; Chakravarthy and Gupta 2020).

	N
	64	128	256	512	1024	2048
$E_1^N$	1.6443e-03	8.2997e-04	4.1700e-04	2.0901e-04	1.0463e-04	5.2349e-05
$E_2^N$	9.0260e-04	4.5166e-04	2.2593e-04	1.1299e-04	5.6502e-05	2.8253e-05
(Chakravarthy and Gupta 2020)
$E_1^N$	3.0136e–03	1.5232e–03	7.6577e–04	3.8392e–04	1.9222e–04	9.6175e–05
$E_2^N$	1.8168e–03	9.3052e–04	4.7082e–04	2.3681e–04	1.1875e–04	5.9465e–05
(Selvi and Ramanujam 2017)
$E_1^N$	5.7577e-03	2.9662e-03	1.4985e-03	7.5187e-04	3.7640e-04	1.8828e-04
$E_2^N$	4.5703e-03	2.5914e-03	1.4979e-03	8.3942e-04	4.6874e-04	2.5666e-04
(Subburayan and Ramanujam 2014)
$E_1^N$	3.8838e-03	1.9843e-03	9.9963e-04	5.0116e-04	2.5084e-04	1.2548e-04
$E_2^N$	9.8093e-03	6.9348e-03	4.4954e-03	2.8118e-03	1.6521e-03	9.4091e-04

Table 2. Rate of convergence $\mathit{\rho }$ of the example 1 for $\mathit{ϵ}=2^{-8}$.

	$\mathit{h}$	$\frac{h}{2}$	$Z_h$	$\frac{h}{4}$	$Z_{\frac{h}{2}}$	$\mathit{\rho }$
	1/100	1/200	1.0595e-03	1/400	5.3390e-04	0.9888
$R_1^N$	1/200	1/400	5.3390e-04	1/800	2.8150e-04	0.9234
	1/300	1/600	3.6209e-04	1/1200	2.0707e-04	0.8062
	1/100	1/200	5.7799e-04	1/400	2.8923e-04	0.9988
$R_2^N$	1/200	1/400	2.8923e-04	1/800	1.4713e-04	0.9751
	1/300	1/600	1.9357e-04	1/1200	1.0399e-04	0.8965

The graph of the computed solution of two components for different values of $\mathit{\epsilon }$ is given in Figure 1. Figure 3 represents the graphs of maximum absolute error of both components for different values $\mathit{\epsilon }$, whereas Figure 2 represents the graphs of point-wise absolute errors for different values of $\mathit{N}$ with $\mathit{\epsilon }=0.1$. The loglog plot of the maximum pointwise error is also given in Figure 4.

Figure 1. The numerical solution of example 1 with $\mathit{\epsilon }=0.1$ and $\mathit{N}=100$.

Figure 2. The point-wise absolute errors of example 1 for different values of $\mathit{N}$, with $\mathit{\epsilon }=0.1$.

Figure 3. The maximum absolute error of example 1 for different values of $\mathit{\epsilon }$.

Figure 4. Loglog plot of maximum point-wise errors of example 1 for different values of $\mathit{\epsilon }$.

6.2.

Example 2

$\left\{\begin{array}{ccc}-\mathit{\epsilon }{\mathit{z}{\mathit{ }}^{\mathrm{\prime \prime }}}_1(\mathit{t})+11z_{1{\mathit{ }}^{\mathrm{\prime }}}(\mathit{t})+6z_1(\mathit{t})-2z_2(\mathit{t})-z_1(\mathit{t}-1)& =& 0,\mathit{t}\in [0,2],\\ -\mathit{\epsilon }{\mathit{z}{\mathit{ }}^{\mathrm{\prime \prime }}}_2(\mathit{t})+16z_{2{\mathit{ }}^{\mathrm{\prime }}}(\mathit{t})-2z_1(\mathit{t})+5z_2(\mathit{t})-z_2(\mathit{t}-1)& =& 0,\mathit{t}\in [0,2],\\ z_1(\mathit{t})=1,\mathit{t}\in [-1,0],z_1(2)=1,& \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}\\ z_2(\mathit{t})=1,\mathit{t}\in [-1,0],z_2(2)=1.& \mathit{\hspace{0.17em}}& \mathit{\hspace{0.17em}}\end{array}\right.$

The maximum pointwise errors are presented in Table 4 for different values of $\mathit{\epsilon }$ and the numerical results are compared with the results in (Selvi and Ramanujam 2017). As the step size $\mathit{h}=1/\mathit{N}$ tends to zero for all values of the perturbation parameter epsilon, the maximum absolute errors decrease with the linear rate of convergence, as shown by the results presented in Table 5 for different values of $\mathit{N}$ with $\mathit{\epsilon }=2^{-8}$.

Table 4. The maximum absolute error of the example 2 and result in (Selvi and Ramanujam 2017).

	N
	64	128	256	512	1024	2048
$E_1^N$	1.1705e-03	6.0307e-04	3.0609e-04	1.5420e-04	7.7390e-05	3.8768e-05
$E_2^N$	8.3655e-04	4.2276e-04	2.1250e-04	1.0653e-04	5.3336e-05	2.6685e-05
(Selvi and Ramanujam 2017)
$E_1^N$	1.0462e-02	6.5658e-03	4.2022e-03	2.5157e-03	1.4694e-03	8.3463e-04
$E_2^N$	7.2444e-03	5.3271e-03	3.5089e-03	2.1973e-03	1.2983e-03	7.4184e-04

Table 5. Rate of convergence $\mathit{\rho }$ of the example 2 for $\mathit{\epsilon }=2^{-8}$.

	$\mathit{h}$	$\frac{h}{2}$	$Z_h$	$\frac{h}{4}$	$Z_{\frac{h}{2}}$	$\mathit{\rho }$
	1/100	1/200	1.0595e-03	1/400	5.3390e-04	0.9888
$R_1^N$	1/200	1/400	5.3390e-04	1/800	2.8150e-04	0.9234
	1/300	1/600	3.6209e-04	1/1200	2.0707e-04	0.8062
	1/100	1/200	5.7799e-04	1/400	2.8923e-04	0.9988
$R_2^N$	1/200	1/400	2.8923e-04	1/800	1.4713e-04	0.9751
	1/300	1/600	1.9357e-04	1/1200	1.0399e-04	0.8965

The graph of the computed solution of two components for different values of $\mathit{\epsilon }$ is given in Figure 5. Figure 6 represents the graphs of maximum absolute error of both components for different values $\mathit{\epsilon }$, whereas Figure 7 represents the graphs of point-wise absolute errors for different values of $\mathit{N}$ with $\mathit{\epsilon }=0.1$. The loglog plot of the maximum pointwise error is also given in Figure 8.

Figure 5. The numerical solution of example 2 with $\mathit{\epsilon }=0.1$ and $\mathit{N}=100$.

Figure 6. The maximum absolute error of example 2 for different values of $\mathit{\epsilon }$.

Figure 7. The point-wise absolute errors of example 2 for different values of $\mathit{N}$, with $\mathit{\epsilon }=0.1$.

Figure 8. Loglog plot of maximum point-wise errors of example 1 for different values of $\mathit{\epsilon }$.

The Rate of Convergence $(R_i^N)$

The computational rate of convergence $R_i^N$ is also obtained by using the double mesh principle defined as:

$R_i^N=\frac{log(E_i^N)-log(E_i^{2N})}{log2},\mathit{i}=1,2$

7. Conclusion

This work presents an effective numerical approach for handling a system of SPDDEs (convection-diffusion problems) with large delay. To demonstrate the efficacy of this technique, various values of $\mathit{\epsilon }$ have been applied to some model examples (which do not have an actual solution). It is apparent from the tables that the methodology provides an approximation to the solution. The numerical solutions are summarized in terms of maximum absolute errors Tables 1 and 4, and it is discovered that the discussed method gives an improvement compared to (Subburayan and Ramanujam 2014; Selvi and Ramanujam 2017; Chakravarthy and Gupta 2020). The findings demonstrate convergence of the obtained solution by showing that the maximum absolute errors decline with decreasing grid size and the convergence is of order one with this method. The graphs represent the numerical solution of the problem. Tables 2 and 5 represent the rate of convergence $R_i^N$.

Acknowledgments

The authors express their sincere thanks to the editor and reviewers for their valuable comments and suggestions that improved the quality of the manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

No data was used for the research described in the article.

Additional information

Funding

This work was supported by the VIT University.