Numerical Treatment of the Coupled Fredholm Integro-Differential Equations by Compact Finite Difference Method

Numerical Treatment of the Coupled Fredholm Integro-Differential Equations by Compact Finite Difference Method

Surme R. Saber Younis A. Sabawi* Hoshman Q. Hamad Mohammad Sh. Hasso

Department of Mathematics, Faculty of Science and Health, Koya University, Koya 46017, Iraq

College of Information Technology, Imam Ja’afar AL-Sadiq University, Baghdad 10011, Iraq

Corresponding Author Email: 
younis.abid@koyauniversity.org
Page: 
2671-2685
|
DOI: 
https://doi.org/10.18280/mmep.111009
Received: 
29 November 2023
|
Revised: 
28 April 2024
|
Accepted: 
10 May 2024
|
Available online: 
31 October 2024
| Citation

© 2024 The authors. This article is published by IIETA and is licensed under the CC BY 4.0 license (http://creativecommons.org/licenses/by/4.0/).

OPEN ACCESS

Abstract: 

The work introduces a novel numerical method for solving the Fredholm Integro-Differential Equations (FIDEs) and system of Fredholm Integro-Differential Equations (SFIDEs) by employing the fourth-order compact finite difference methods in conjunction with Simpson’s quadrature rule. The accuracy of the proposed scheme is rigorously evaluated using $l^2$ and $l^{\infty}$ norms, while the computational efficiency is measured by assessing the CPU-time values, demonstrating a notable reduction in computational cost compared to standard finite difference schemes. The significance of this approach lies in its ability to maintain high levels of accuracy, addressing a common challenge in traditional methods. The methods presented exhibit fourth accuracy in space, as evidenced by numerical experiments. The mentioned work signifies a notable progress in tackling problems related to FIDEs and SFIDEs. It introduces a robust and efficient numerical methodology that proves particularly effective in situations where obtaining exact solutions poses challenges. This advancement is crucial as it addresses a common difficulty faced in the solution of FIDEs and SFIDEs problems, offering a reliable numerical approach that can handle complex scenarios and contribute to more accurate and practical solutions in various fields of study.

Keywords: 

Fredholm Integro-Differential Equations, Compact Finite Difference Method, quadrature method

1. Introduction

Recent years, interest is particularly pronounced in the domains of physical and biological modeling, as well as in (bio-)engineering applications where initial- and boundary-value problems involving Integro-Differential Equations are prevalent. The Compact Finite Difference Method (CFDM) has emerged as a valuable tool for approximating solutions to these intricate problems. Its application extends to both ordinary and partial differential equations, making it a versatile and sought-after numerical technique. The compact nature of the method, coupled with its ability to handle Integro-Differential Equations, contributes to its growing appeal in various scientific and engineering disciplines. Researchers and practitioners are increasingly turning to CFDM as an effective approach to address the challenges posed by complex systems described by Integro-Differential Equations, fostering advancements in computational methodologies and facilitating more accurate modelling in diverse fields [1-9].

In light of the limited progress in the development of high-order compact finite difference methods for Integro-Differential Equations (IDE), this research is dedicated to addressing this gap by focusing on the derivation of efficient and accurate high-order numerical schemes. The complexity inherent in Integro-Differential Equations, particularly in the context of second-order SFIDEs, necessitates advanced numerical tools for precise approximations of solutions. By concentrating on second-order SFIDEs, this work aims to contribute significantly to the ongoing efforts in developing sophisticated numerical methodologies tailored for Integro-Differential Equations. The research endeavours to provide a comprehensive understanding of the challenges posed by second-order SFIDEs and to propose novel high-order compact finite difference methods as a solution. The outcomes of this study have the potential to not only advance the field of numerical analysis but also to enhance the capabilities for tackling real-world problems in diverse scientific and engineering applications [10-12].

Consider the second-order FIDEs as follows:

$\begin{gathered}u^{\prime \prime}(x)+q(x) u^{\prime}(x)+p(x) u(x)=f(x)+\lambda \int_a^b k(x, t) u(t) d t\end{gathered}$          (1)

with Dirichlet boundary conditions:

$u(a)=\alpha, \quad u(b)=\beta$          (2)

A system of second-order SFIDEs can be defined as:

$\begin{gathered}u^{\prime \prime}(x)+q_1(x) v^{\prime}(x)+p_1(x) u(x)=f_1(x)+\lambda_1 \int_a^b\left(k_1(x, t) u(t)+k_2(x, t) v(t)\right) d t\end{gathered}$          (3)

$\begin{gathered}v^{\prime \prime}(x)+q_2(x) u^{\prime}(x)+p_2(x) v(x)=f_2(x)+\lambda_2 \int_a^b\left(k_3(x, t) u(t)+k_4(x, t) v(t)\right) d t\end{gathered}$,          (4)

with Dirichlet boundary conditions:

$u(a)=\alpha_1,\ u(b)=\beta_1,\ v(a)=\alpha_2,\ v(b)=\beta_2$.          (5)

For $x,\ t \in[a, b]$, where all $\lambda, \alpha$, and $\beta$ are constant values, all functions $f(x),\ q(x),\ p(x)$ and $k(x, t)$ are known and $u(x)$ and $v(x)$ are the solutions to be determined.

In recent times, there has been a notable surge in interest regarding the application of higher-order numerical methods to integro-differential equations. A particularly promising approach gaining attention is the utilization of compact difference approximations, known for their ability to achieve high accuracy even with a relatively small number of grid points. This method strategically employs 9 grid points, forming a compact patch of three cells around a chosen node, to effectively nullify second-order truncation error terms. The key advantage lies in the subsequent development of alternative expressions, with lower derivatives, that are equivalent to the higher-order truncation error terms. This innovative strategy holds the potential to yield more efficient and accurate solutions for partial differential equations (PDEs), presenting a compelling avenue for advancing the field of numerical analysis and enhancing the precision of solutions in various scientific and engineering applications.

It is known that compact difference approximations exist for certain operators that are higher-order than standard schemes. As an example, in the reference [13], the difference coefficients at the five grid points corresponding to the compact patch of three cells surrounding a given node can be selected so the second-order truncation error terms cancel. Therefore, utilizing the differential equation, leads to producing alternative lower-derivative expressions equivalent to the higher-order truncation error terms [14].

The primary objective of this research is to establish a comprehensive framework and methodology for the development of higher-order compact (HOC) schemes tailored for second-order System of Fredholm Integro-Differential Equations (SFIDEs). Drawing inspiration from various studies employing compact finite difference methods for Volterra integral problems, including references [15-17]. Numerous authors have developed numerical methods for integral and Integro-Differential Equations recently, see references [18-24].

The present study is prepared as follows: In Section 2, the numerical method for second-order FIDEs and the system of second-order FIDEs equations are proposed. In Section 3, convergence analysis is proposed. Some numerical experiments are shown in Section 4. Finally, conclusions are given in Section 5.

2. Fourth-Order Compact Finite Difference Scheme

2.1 Second-order Fredholm Integro-Differential Equations

We employ the fourth-order accuracy Compact Finite Difference Method to tackle the problem stated in Eq. (1). To begin, discretize the interval $I=[a, b]$ into N equal spaced grid points with nodes $x_0, x_1, x_2, \ldots x_{N-1}, x_N$, then

$I=\left\{a=x_0<x_1<x_2<\ldots<x_{N-1}<x_N=b\right\}$,

and $x_i=x_0+i h, i=0,1,2, \ldots, N$. ($h$ is the spatial step size). The following notation $u(x+n h)=u_{i+n}$ at $x=x_i$ is used for clarity.

By Taylor’s series expansion

$\begin{aligned} u_{i+1} & =u_i+h u_i^{\prime}+\frac{h}{2!} u_i^{\prime \prime}+\frac{h^3}{3!} u_i^{(3)}+\frac{h^4}{4!} u_i^{(4)}+\frac{h^5}{5!} u_i^{(5)}+\frac{h^6}{6!} u_i^{(6)}+\frac{h^7}{7!} u_i^{(7)}+\cdots\end{aligned}$         (6)

$u_{i-1}=u_i-h u_i^{\prime}+\frac{h}{2!} u_i^{\prime \prime}-\frac{h^3}{3!} u_i^{(3)}+\frac{h^4}{4!} u_i^{(4)}-\frac{h^5}{5!} u_i^{(5)}+\frac{h^6}{6!} u_i^{(6)}-\frac{h^7}{7!} u_i^{(7)}+\cdots$.           (7)

Denoting the central difference of second order approximate with first and second derivative of the function $u(x)$ by $\delta_x u_i$ and $\delta_x^2 u_i$ respectively, where $\delta_x u_i$ and $\delta_x^2 u_i$ are given as

$\delta_x u_i=\frac{u_{i+1}-u_{i-1}}{2 h}, \delta_x^2 u_i=\frac{u_{i+1}-2 u_i+u_{i-1}}{h^2}$.           (8)

Substituting Taylor's series expansion for $u_{i+1}$ and $u_{i-1}$ in (6) and (7) into (8) we have that

$\begin{aligned} \delta_x u_i & =u_i^{\prime}+\frac{h^2}{3!} u_i^{(3)}+\frac{h^4}{5!} u_i^{(5)}+\frac{h^6}{7!} u_i^{(7)}+\frac{h^8}{9!} u_i^{(9)}+\frac{h^{10}}{11!} u_i^{(11)}+\cdots\end{aligned}$           (9)

$\begin{aligned} \delta_x^2 u_i=u_i^{\prime \prime}+ & 2\left(\frac{h^2}{4!} u_i^{(4)}+\frac{h^4}{6!} u_i^{(6)}+\frac{h^6}{8!} u_i^{(8)}\right.\left.+\frac{h^8}{10!} u_i^{(10)}+\frac{h^{10}}{12!} u_i^{(12)}+\cdots\right)\end{aligned}$.           (10)

So that

$\begin{aligned} u_i^{\prime}=\delta_x u_i- & \left(\frac{h^2}{3!} u_i^{(3)}+\frac{h^4}{5!} u_i^{(5)}+\frac{h^6}{7!} u_i^{(7)}\right.\left.+\frac{h^8}{9!} u_i^{(9)}+\frac{h^{10}}{11!} u_i^{(11)}+\cdots\right)\end{aligned}$           (11)

$\begin{aligned} u_i^{\prime \prime}=\delta_x^2 u_i- & 2\left(\frac{h^2}{4!} u_i^{(4)}+\frac{h^4}{6!} u_i^{(6)}+\frac{h^6}{8!} u_i^{(8)}\right.\left.+\frac{h^8}{10!} u_i^{(10)}+\frac{h^{10}}{12!} u_i^{(12)} \cdots\right)\end{aligned}$,           (12)

substituting (11) and (12) into problem (1) yields

$L_h u_i-\operatorname{Err} 1=f_i+\lambda \int_a^b k_i(t) u(t) d t$           (13)

where,

$L_h u_i=\delta_x^2 u_i+q_i \delta_x u_i+p_i u_i$,

where, Err1 is the truncation error given as

$\begin{aligned}

& \operatorname{Err} 1=2\left(\frac{h^2}{4!} u_i^{(4)}+\frac{h^4}{6!} u_i^{(6)}+\frac{h^6}{8!} u_i^{(8)}\right)+\left(\frac{h^2}{3!} u_i^{(3)}+\frac{h^4}{5!} u_i^{(5)}+\frac{h^6}{7!} u_i^{(7)}\right)+\cdots \\

& \operatorname{Err} 1=\frac{2 h^2\left(2 q_i u_i^{(3)}+u_i^{(4)}\right)}{4!} +\frac{2 h^4\left(3 q_i u_i^{(5)}+u_i^{(6)}\right)}{6!}+\frac{2 h^6\left(4 q_i u_i^{(7)}+u_i^{(8)}\right)}{8!}+\cdots,

\end{aligned}$               (14)

where, Err1=$E_4+E_6+E_8+\cdots$.

To achieve 4th-order accuracy, we can focus on the error term $E_4$ in the expansion of the central difference of the second order approximation given in (4). By summing this error term to the central difference of the second order scheme, we can obtain a scheme with 4th-order accuracy such that

$E_4=\frac{h^2}{12}\left[2 q_i u_i^{(3)}+u_i^{(4)}\right]$.           (15)

To obtain the terms $u_i^{(3)}$ and $u_i^{(4)}$, we begin by rearranging Eq. (1) and then differentiating it repeatedly with respect to x from Eq. (1), we have

$\boldsymbol{u}_i^{\prime \prime}=-q_i u_i^{\prime}-p_i u_i+f_i+\lambda \int_a^b k_i(t) u(t) d t$,             (16)

$\boldsymbol{u}_i^{(3)}=-q_i^{\prime} u_i^{\prime}-q_i u_i^{\prime \prime}-p_i^{\prime} u_i-p_i u_i^{\prime}+f_i^{\prime}+\int_a^b k_i^{\prime} u(t) d t$               (17)

$\boldsymbol{u}_i^{(4)}=\left(q_i^2-2 q_i^{\prime}-p_i\right) u_i^{\prime \prime}+\left(q_i q_i^{\prime}+q_i p_i-q_i^{\prime \prime}-2 p_i\right) u_i^{\prime}+\left(q_i p_i^{\prime}-p_i^{\prime \prime}\right) u_i-q_i f_i^{\prime}+f_i^{\prime \prime}-q_i \int_a^b k_i^{\prime}(t) t u(t) d t+\int_a^b k_i^{\prime \prime}(t) u(t) d t$.             (18)

Then substituting $\boldsymbol{u}_i^{(3)}$ from (17) and $\boldsymbol{u}_i^{(4)}$ from (18) into the term $E_4$ in (15) and simplifying produces

$\begin{aligned} E_4 & =-\left(\frac{h^2\left(q_i^2+2 q_i^{\prime}+p_i\right)}{12}\right) u_i^{\prime \prime}-\left(\frac{h^2\left(q_i q_i^{\prime}+q_i p_i+q_i^{\prime \prime}+2 p_i^{\prime}\right)}{12}\right) u_i^{\prime}-\left(\frac{h^2\left(q_i p_i^{\prime}-p_i^{\prime \prime}\right)}{12}\right) u_i \\ & +\left(\frac{h^2 q_i}{12}\right) f_i^{\prime}+\frac{h^2}{12} f_i^{\prime \prime}+\left(\frac{h^2 q_i}{12}\right) \int_a^b k_i^{\prime}(t) u(t) d t+\left(\frac{h^2}{12}\right) \int_a^b k_i^{\prime \prime}(t) u(t) d t .\end{aligned}$.            (19)

Replace the terms $u_i^{\prime \prime}$ and $u_i^{\prime}$ in $E_4$ with those in (11) and (12) to obtain

$\begin{aligned} E_4 & =-\left(\frac{h^2\left(q_i^2+2 q_i^{\prime}+p_i\right)}{12}\right) \delta_x^2 u_i-\left(\frac{h^2\left(q_i q_i^{\prime}+q_i p_i+q_i^{\prime \prime}+2 p_i^{\prime}\right)}{12}\right) \delta_x u_i-\left(\frac{h^2\left(q_i p_i^{\prime}-p_i^{\prime \prime}\right)}{12}\right) u_i \\ & +\left(\frac{h^2 q_i}{12}\right) f_i^{\prime}+\frac{h^2}{12} f_i^{\prime \prime}+\left(\frac{h^2 q_i}{12}\right) \int_a^b k_i^{\prime}(t) u(t) d t+\left(\frac{h^2}{12}\right) \int_a^b k_i^{\prime \prime}(t) u(t) d t-E r r 2,\end{aligned}$       (20)

the error term Err2 is given as

$\begin{gathered}\text { Err } 2=2\left(\frac{h^2\left(q_i^2+2 q_i^{\prime}+p_i\right)}{12}\right) \\ \left(\frac{h^2}{4!} u_i^{(4)}+\frac{h^4}{6!} u_i^{(6)}+\frac{h^6}{8!} u_i^{(8)}+\frac{h^8}{10!} u_i^{(10)}+\frac{h^{10}}{12!} u_i^{(12)} \cdots\right)+\left(\frac{h^2\left(q_i q_i^{\prime}+q_i p_i+q_i^{\prime \prime}+2 p_i^{\prime}\right)}{12}\right) \\ \left(\frac{h^2}{3!} u_i^{(3)}+\frac{h^4}{5!} u_i^{(5)}+\frac{h^6}{7!} u_i^{(7)}+\frac{h^8}{9!} u_i^{(9)}+\frac{h^{10}}{11!} u_i^{(11)}+\cdots\right)\end{gathered}$

Adding $E_4$ to the second order central difference scheme, gives

$\begin{gathered}\left(1+\frac{h^2\left(q_i^2+2 q_i^{\prime}+p_i\right)}{12}\right) \delta_x^2 u_i+\left(q_i+\frac{h^2\left(q_i q_i^{\prime}+q_i p_i+q_i^{\prime \prime}+2 p_i^{\prime}\right)}{12}\right) \delta_x u_i+\left(\frac{h^2\left(q_i p_i^{\prime}-p_i^{\prime \prime}\right)}{12}\right) u_i+\left(\frac{h^2 q_i}{12}\right) f_i^{\prime}+\frac{h^2}{12} f_i^{\prime \prime} \\ +\left(\frac{h^2 q_i}{12}\right) \int_a^b k_i^{\prime}(t) u(t) d t+\left(\frac{h^2}{12}\right) \int_a^b k_i^{\prime \prime}(t) u(t) d t-{Err} 3=f_i+\lambda \int_a^b k_i(t) u(t) d t .\end{gathered}$       (21)

The new error term Err3 is a linear combination of the error obtained from the central difference scheme (8) and the error obtained from $E_4$ i.e. Err3=Err1+Err2, so that

$\begin{aligned} \operatorname{Err} 3 & =h^4\left(\frac{2\left(3 q_i u_i^{(5)}+u_i^{(6)}\right)}{6!}+\frac{k_1}{3!} u_i^{(3)}+\frac{2 k_2}{4!} u_i^{(4)}\right)+h^6\left(\frac{2\left(4 q_i u_i^{(7)}+u_i^{(8)}\right)}{8!}+\frac{k_1}{5!} u_i^{(5)}+\frac{2 k_2}{6!} u_i^{(6)}\right) \\ & +h^8\left(\frac{2\left(3 q_i u_i^{(9)}+u_i^{(10)}\right)}{10!}+\frac{k_1}{7!} u_i^{(7)}+\frac{2 k_2}{8!} u_i^{(8)}\right)+h^{10}\left(\left(\frac{2\left(3 q_i u_i^{(11)}+u_i^{(12)}\right)}{12!}+\frac{k_1}{9!} u_i^{(9)}+\frac{2 k_2}{10!} u_i^{(10)}\right)\right)+\cdots ,\end{aligned}$

where, $k_1=\frac{q_i q_i^{\prime}+q_i p_i+q_i^{\prime \prime}+2 p_i^{\prime}}{12}$ and $k_1=\frac{q_i^2+2 q_i^{\prime}+p_i}{12}$.

The integral part on the right-hand side of Eq. (21) will be handled numerically using the composite Simpson’s rule given by:

$I_i=\Upsilon_i^h(k u)+\mathcal{T}^{S M}$,

where,

$\Upsilon_i^h(k u)=\frac{h}{3}\left(k_{i, 0} u_0+4 \sum_{j=1}^{\frac{n}{2}} k_{i, 2 j-1} u_{2 j-1}+2 \sum_{j=1}^{\left(\frac{n}{2}\right)-1} k_{i, 2 j} u_{2 j}+k_{i, n} u_n\right)$,          (22)

and $\mathcal{T}^{S M}$ be a vector denoting the truncation error such that where $\left\|\mathcal{T}^{S M}\right\|=O\left(h^4\right)$. Replacing the terms $\delta_x u_i$ and $\delta_x^2 u_i$ in (22) with their definitions in (4) and neglecting the error term produces following fourth-order compact finite difference scheme

$\begin{gathered}A_i u_{i+1}+B_i u_i+C_i u_{i-1}=F_i+I_i \\ i=1,2,3, \ldots, N-1,\end{gathered}$           (23)

where, the coefficient $A_i, B_i, C_i, F_i$ and $I_i$ are given as

$A_i=\frac{1}{h^2}+\frac{\left(q_i^2+2 q_i^{\prime}+p_i\right)}{12}+\frac{q_i}{2 h}+\frac{h\left(q_i q_i^{\prime}+q_i p_i+q_i^{\prime \prime}+2 p_i^{\prime}\right)}{24}$,

$B_i=\frac{-2}{h^2}+\frac{h\left(q_i^2+2 q_i^{\prime}+p_i\right)}{6}+\frac{h^2\left(q_i p_i^{\prime}-p_i^{\prime \prime}\right)}{12}$,

$C_i=\frac{1}{h^2}+\frac{\left(q_i^2+2 q_i^{\prime}+p_i\right)}{12}-\frac{q_i}{2 h}-\frac{h\left(q_i q_i^{\prime}+q_i p_i+q_i^{\prime \prime}+2 p_i^{\prime}\right)}{24}$,

$F_i=f_i-\left(\frac{h^2 q_i}{12}\right) f_i^{\prime}-\frac{h^2}{12} f_i^{\prime \prime}$,

$I_i=\lambda \int_a^b k_i(t) u(t) d t+\left(\frac{h^2 q_i}{12}\right) \int_a^b k_i^{\prime} u(t) d t-\left(\frac{h^2}{12}\right) \int_a^b k_i^{\prime \prime} u(t) d t$,

where, $u_0$ and $u_n$ are Dirichlet boundary conditions get in general system, when i=1, and i=n-1 we obtain $u_0$ and $u_n$ respectively. The system in Eq. (23) consists of (N-1) linear equations with (N-1) unknowns $\left(u_1, u_2, u_3, \ldots, u_{n-2}, u_{n-1}\right)$.

2.2 The algorithm

In order to find the numerical solution of Eq. (1) by utilizing the fourth-order compact difference approximation with the composite Simpson’s rule, we provide the method of solution in algorithm

Algorithm 1

Input: N, a, b and boundary condition $u_0$ and $u_n$, where $\left(a=x_0, b=x_n\right)$.

Set: $h=\frac{b-a}{N}$.

for $i \leftarrow 0$ to N do

     for $j \leftarrow 0$ to N do

           $x_i=a+i h$.

           $t_j=a+j h$.

     end for

end for

for $i \leftarrow 0$ to N do

           $u_i \leftarrow u\left(x_i\right)$

          for $j \leftarrow 0$ to N do

               $k_{i, j} \leftarrow k\left(x_i, t_j\right)$

          end for

end for

for $i \leftarrow 1$ to N-1 do

     $C_i=\mathcal{F}+\omega$.

     for $j \leftarrow 1$ to N-1 do

          $\mathcal{A}_{(i, j)}=\mathcal{B}+\gamma$.

     end for

end for

$C=\left[C_1-\gamma_1 u_0 ; C_2: C_{n-2} ; C_{n-1}-\gamma_3 u_n\right]$

Output: $U \leftarrow \mathcal{A} \backslash C$

3. Convergence Analysis

This section aims to prove the convergence analysis of the fourth-order CFD method. We first introduce the following lemma that play a crucial role in the proof. To prove the convergence of the system in Eq. (36) using the same method as the convergence analysis of the fourth-order CFD method. We define the space $L_2(a, b)$ represents a Hilbert space with the inner product

$(u(x), v(x))=\int_a^b u(x) . v(x) d x$.

The Sobolev $H^1$ admit a natural norm

$\|u\|_{H^1}^2=\|u\|^2+\left\|u^{\prime}\right\|^2$.

Lemma 3.1: The remainder $\mathcal{T}^{S M}$ of the composite Simpson’s rule satisfies

$\left|\mathcal{T}^{S M}\right| \leq \frac{1}{180} h^4 u^{(i v)}(\xi)$.

Proof. Let the function $z=K\left(x_i, s\right) u(s)$ be continuous and possess a continuous derivative in $\left[s_0, s_2\right]$. Expanding y about $s=s_0$ we obtain

$\begin{gathered}z(x)=z_0+\left(s-s_0\right) z_0^{\prime}+\frac{1}{2}\left(s-s_0\right)^2 z_0^{\prime \prime}+\frac{1}{3!}\left(s-s_0\right)^3 z_0^{\prime \prime \prime}+\frac{1}{4!}\left(s-s_0\right)^3 z_0^{(i v)}+\cdots \int_{s_0}^{s_2} K\left(x_i, s\right) u(s) d s \\ =\int_0^2 h\left(z_0+r h z_0^{\prime}+\frac{(r h)^2}{2} z_0^{\prime \prime}+\frac{(r h)^3}{3!} z_0^{\prime \prime \prime}+\frac{(r h)^4}{4!} z_0^{(i v)}+\cdots\right) d r \\ =h\left[r z_0+\frac{r^2 h}{2} z_0^{\prime}+\frac{(r h)^2}{6} p z_0^{\prime \prime}+\frac{(r h)^3}{24} p z_0^{\prime \prime \prime}+\frac{(r h)^4}{120} r z_0^{(i v)}+\right]_0^2 \\ =2 h z_0+2 h^2 z_0^{\prime}+\frac{4 h^2}{3} z_0^{\prime \prime}+\frac{2 h^4}{3} z_0^{\prime \prime \prime}+\frac{4 h^5}{15} z_0^{(i v)}+\cdots.\end{gathered}$          (24)

Therefore,

$z_0=z_0$        (25)

$z_1=z_0+h z_0^{\prime}+\frac{h^2}{2} z_0^{\prime \prime}+\frac{h^3}{6} z_0^{\prime \prime \prime}+\frac{h^4}{24} z_0^{(i v)}+\cdots$,         (26)

$z_2=z_0+2 h z_0^{\prime}+2 h^2 z_0^{\prime \prime}+\frac{4 h^3}{3} z_0^{\prime \prime \prime}+\frac{2 h^4}{3} z_0^{(i v)}+\cdots$.          (27)

Combining (24), (25) and (26), becomes

$\frac{h}{3}\left[z_0+4 z_1+z_2\right]=\frac{h}{3}\left[6 z_0+6 h z_0^{\prime}+4 h^2 z_0^{\prime \prime}+2 h^3 z_0^{\prime \prime \prime}+\frac{5 h^4}{6} z_0^{(i v)}+\cdots,\right]=2 h z_0+2 h^2 z_0^{\prime}+\frac{4 h^3}{3} z_0^{\prime \prime}+\frac{2 h^4}{3} z_0^{\prime \prime \prime}+\frac{5 h^4}{18} z_0^{(i v)}+\cdots$.              (28)

Using (24) and (28), this leads

$\int_{s_0}^{s_2} z d s-\frac{h}{3}\left[z_0+4 z_1+z_2\right]=\frac{-1}{90} h^5 z_0^{(i v)}$.

The composite Simpson’s 1/3 rule is a numerical integration method used to approximate the definite integral of a function over an interval by dividing the interval into multiple sub-intervals and applying Simpson’s 1/3 rule on each sub-interval. To use this method, the interval of integration [0,l] is subdivided into N even number of sub-divisions as follows: $0=s_0<s_1<s_2<\cdots<s_N=l$. The integral over the entire interval [0,l] can be approximated by summing up the individual integrals over each sub-interval using the composite Simpson’s 1/3 rule. The formula for the composite Simpson’s 1/3 rule for a sub-interval $\left[s_i, s_{i+1}\right]$ is given by:

$\int_{s_0}^{s_2} z d s=\frac{h}{3}\left(z_0+4 \sum_{j=1}^{N / 2} z_{2 j-1}+2 \sum_{j=1}^{(N / 2)-1} z_{2 j}+z_N.\right)=\frac{-1}{90} h^5 z_0^{(i v)}$.

We obtain the errors in the intervals [0,l] as

$\mathcal{T}^{S M}=\frac{-1}{90} h^5\left[z_0{ }^{(i v)}+z_2{ }^{(i v)}+\cdots+z_{N-2}{ }^{(i v)}\right]=\frac{-1}{180} h^4 u^{(i v)}(\xi), \xi \in[0, l]$.

Assumption 1: The kernel functions k(x,s) satisfies the following positive definite property:

$\int_a^b \int_a^b(k(x, s) \theta(x), \theta(s)) d x d s>0$.           (29)

for every continuous $\theta(x)=\left(\theta_1(x), \theta_2(x), \ldots \theta_k(x)\right) \neq 0$, and the integral

$\int_a^b \int_a^b|k(x, s)|^2 d x \,\, d s<\infty$.

Theorem 3.2 (Error estimates). Assuming that assumption 1 is satisfied and the kernel function k(x, s) is smooth enough, then

$\left\|u(x)-\Theta_h(x)\right\|_{H^1} \leq C h^4$,          (30)

where $C>0$.

Proof: Setting $e_h=u-\Theta_h(x)$ in Eq. (1), gives

$e_h^{\prime \prime}+\mathrm{p}(\mathrm{x}) e_h^{\prime}+\mathrm{q}(\mathrm{x}) e_h-\gamma \int_a^b k(\mathrm{x}, s) e_h(s) d s=\mathrm{f}(x)+\mathcal{T}(x, h)$,           (31)

where $\mathcal{T}(x, h)$ be a vector denoting the truncation error such that

$\mathcal{T}(x, h)=-\Theta_h^{\prime \prime}-\mathrm{p}(\mathrm{x}) \Theta_h^{\prime}-\mathrm{q}(\mathrm{x}) \Theta_h+\lambda \int_a^b k(\mathrm{x}, s) \Theta_h(s) d s$.             (32)

Then, we have

$\mathcal{T}\left(x_i, h\right)=-L_h z_i+\lambda \int_a^b k\left(x_i, s\right) \mathrm{w}(z, s) d s-\Theta_h^{\prime \prime}-\mathrm{p}\left(x_i\right) \Theta_h^{\prime}-\mathrm{q}\left(x_i\right) \Theta_h+L_h z_i-\lambda \int_a^b k(x, s)\left(\mathrm{w}(z, s)-\Theta_h(x)\right) d x=\mathrm{O}\left(h^4\right)$.         (33)

Multiplying Eq. (31) by $e_h$ and integrating with respect to x, gives

$\begin{gathered}\int_a^b\left(e_h^{\prime \prime}, e_h\right) d x+\int_a^b q(x)\left(e_h^{\prime}, e_h\right) d x+\int_a^b p(x)\left(e_h, e_h\right) d x \\ -\lambda \int_a^b \int_a^b\left(k(x, s)\left(e_h, e_h\right)\right) d s d x=\int_a^b\left(\mathrm{f}(x), e_h\right) d x+\int_a^b\left(\mathcal{T}(x, h), e_h\right) d x .\end{gathered}$           (34)

Since $p-\frac{1}{2} q^{\prime \prime} \geq 0$ and integrating by part along with Assumption 1, we have

$\int_a^b\left(e_h^{\prime}, e_h^{\prime}\right) d x \leq \int_a^b\left(\mathrm{f}(x), e_h\right) d x+\int_a^b\left(\mathcal{T}(x, h), e_h\right) d x$.          (35)

Applying Cauchy’s inequality and Lemma 3.1 with, truncation error estimate defined in (19), this becomes

$\left\|e_h\right\|_{H^1}^2 \leq\|\mathrm{f}(x)\|\left\|e_h\right\|_{H^1}\|\mathcal{T}(x, h)\|\left\|e_h\right\|_{H^1}$,

which implies

$\left\|e_h\right\|_{H^1} \leq C\|\mathcal{T}(x, h)\|=\left\|\mathcal{T}^{S M}+\operatorname{Err} 4\right\|=C\left\|\mathcal{T}^{S M}+\operatorname{Err} 4\right\| \leq C\left\|\mathcal{T}^{S M}\right\|+\left\|\mathcal{T}^{S M}\right\|\left\|e_h\right\|_{H^1}$

$\leq\left\|\frac{1}{180} h^4 u^{(i v)}(\xi)\right\|+C\left\|h^4\left(\frac{2\left(3 q_i u_i^{(5)}+u_i^{(6)}\right)}{6!}+\frac{k_1}{3!} u_i^{(3)}+\frac{2 k_2}{4!} u_i^{(4)}\right)+h^6\left(\frac{2\left(4 q_i u_i^{(7)}+u_i^{(8)}\right)}{8!}+\frac{k_1}{5!} u_i^{(5)}+\frac{2 k_2}{6!} u_i^{(6)}\right)\right\|$.

$\leq\left\|\frac{1}{180} h^4 u^{(i v)}(\xi)\right\|+\left\|h^4\left(\frac{2\left(3 q_i u_i^{(5)}+u_i^{(6)}\right)}{6!}+\frac{k_1}{3!} u_i^{(3)}+\frac{2 k_2}{4!} u_i^{(4)}\right)\right\|+\left\|h^6\left(\frac{2\left(4 q_i u_i^{(7)}+u_i^{(8)}\right)}{8!}+\frac{k_1}{5!} u_i^{(5)}+\frac{2 k_2}{6!} u_i^{(6)}\right)\right\|$.

By taking small values of $h^4>h^6$, the proof will be finished.

4. Numerical Experiments

In this section, the performance of a method is demonstrated using an implementation in Malab programming. The error norms of $l^2$ and $l^{\infty}$ are utilized to quantify the error between the numerical and analytical solutions. We denote by $E$ errors terms given by:

$E(x)=u(x)-U_{\text {Appro. }}(x)$

  • The pointwise error,

$\varepsilon(x)=\left|E\left(x_i\right)\right|$

  • The $l^{\infty}$-norm of the error,

$l^{\infty}(\mathrm{E}, h)=\max _{0 \leq i \leq N}\left|\mathrm{E}\left(x_i\right)\right|$

  • The $l^2$-norm of the error,

$l^2(\mathrm{E}, h)=\sqrt{h \sum_{i=0}^N\left|\mathrm{E}\left(x_i\right)\right|^2}$

  • The order of convergence R is calculated as,

$Rate=\frac{\log\left(\operatorname{Error}\left(N_1\right) / \operatorname{Error}\left(\mathrm{N}_2\right)\right)}{\log \left(\mathrm{N}_2 / N_1\right)}$

Example 1: We firstly deal with p(x) and q(x)=0 such that:

$u^{\prime \prime}(x)=e^x-x+\int_0^1 x t u(t) d t, 0 \leq x \leq 1$

with Dirichlet boundary conditions: $u(0)=1, u(1)=e^1$, and the exact solution is $u(x)=e^x$.

Example 2: We secondary deal with $p(x)$ and $q(x)=0$ such that:

$u^{\prime \prime}(x)=x-2+60 \int_0^1(x-t) u(t) d t, 0 \leq x \leq 1$

with Dirichlet boundary conditions: $u(0)=0, u(1)=0$, and the exact solution is $u(x)=x(x-1)^2$.

Example 3: We deal with $p(x) \neq 0$ and $q(x)=0$ such that:

$u^{\prime \prime}(x)=9 u(x)+\frac{e^{-15}-1}{3}+\int_0^5 u(t) d t, 0 \leq x \leq 5$

with Dirichlet boundary conditions: $u(0)=0, u(1)=e^{-15}$, and the exact solution is $u(x)=e^{-3 x}$.

4.1 Second-order system Fredholm Integro-Differential Equations

The fourth-order Compact Finite Difference Method, along with the first and second central differentiation formulas defined in Eqs. (6) and (7), provides a high-accuracy approximation for solving differential equations, gives:

$\begin{aligned} & \delta_x^2 u_i+q_{1_i} \delta_x v_i+p_{1_i} u_i-\varepsilon r r 1=f_{1_i}+\lambda_1 \int_a^b\left(k_{1_i}(t) u(t)+k_{2_i}(t) v(t)\right) d t, \\ & \delta_x^2 v_i+q_{2_i} \delta_x u_i+p_{2_i} v_i-\varepsilon r r 2=f_{2_i}+\lambda_2 \int_a^b\left(k_{3_i}(t) u(t)+k_{4_i}(t) v(t)\right) d t,\end{aligned}$         (36)

where, Err1 and Err2 is the truncation error given as

$\begin{gathered}\varepsilon r r 1=2\left(\frac{h^2}{4!} u_i^{(4)}+\frac{h^4}{6!} u_i^{(6)}+\frac{h^6}{8!} u_i^{(8)}\right)+q_{1 i}\left(\frac{h^2}{3!} v_i^{(3)}+\frac{h^4}{5!} v_i^{(5)}+\frac{h^6}{7!} v_i^{(7)}\right)+\cdots \\ \varepsilon r r 2=2\left(\frac{h^2}{4!} v_i^{(4)}+\frac{h^4}{6!} v_i^{(6)}+\frac{h^6}{8!} v_i^{(8)}\right)++q_{2_i}\left(\frac{h^2}{3!} u_i^{(3)}+\frac{h^4}{5!} u_i^{(5)}+\frac{h^6}{7!} u_i^{(7)}\right)+\cdots\end{gathered}$          (37a)

$\begin{aligned} & \varepsilon r r1=\frac{2 h^2\left(2 q_{1_i} v_i^{(3)}+u_i^{(4)}\right)}{4!}+\frac{2 h^4\left(3 q_{1_i} v_i^{(5)}+u_i^{(6)}\right)}{6!}+\frac{2 h^6\left(4 q_{1_i} v_i^{(7)}+u_i^{(8)}\right)}{8!}+\cdots \\ & \varepsilon r r 2=\frac{2 h^2\left(2 q_{2_i} u_i^{(3)}+v_i^{(4)}\right)}{4!}+\frac{2 h^4\left(3 q_{2_i} u_i^{(5)}+v_i^{(6)}\right)}{6!}+\frac{2 h^6\left(4 q_{2_i} u_i^{(7)}+v_i^{(8)}\right)}{8!}+\cdots\end{aligned}$          (37b)

the error term Err1 and Err2 can be written as $\varepsilon r r 1=E_{1_4}+E_{1_6}+E_{1_8}+\cdots$ and $\varepsilon r r 2=E_{2_4}+E_{2_6}+E_{2_8}+\cdots$.

To yield 4th-order accuracy, there is need to only sum the term E4 from Err1 to the central difference of second order scheme given in (8) i.e., from (10)

$\begin{aligned} & E_{u_4}=\frac{h^2}{12}\left[2 q_{1_i} v_i^{(3)}+u_i^{(4)}\right] \\ & E_{v_4}=\frac{h^2}{12}\left[2 q_{2_i} u_i^{(3)}+v_i^{(4)}\right]\end{aligned}$           (38)

The terms $u_i^{(3)}, u_i^{(4)}, v_i^{(3)}$ and $v_i^{(4)}$ are obtained by rearranging the Eq. (3) and differentiating repeatedly. From Eq. (3) we have

$\begin{aligned} & \delta_x^2 u_i=-q_{1 i} \delta_x v_i-p_{1_i} u_i+f_{1_i}+\lambda_1 \int_a^b\left(k_{1_i}(t) u(t)+k_{2_i}(t) v(t)\right) d t, \\ & \delta_x^2 v_i=-q_{2_i} \delta_x u_i-p_{2_i} v_i+f_{2_i}+\lambda_2 \int_a^b\left(k_{3_i}(t) u(t)+k_{4_i}(t) v(t)\right) d t .\end{aligned}$          (39)

From whence we get

$\begin{gathered}\boldsymbol{u}_{\boldsymbol{i}}^{(\mathbf{3})}=-q_{1_i}^{\prime} v_i^{\prime}-q_{1_i} v_i^{\prime \prime}-p_{1 i}^{\prime} u_i-p_{1_i} u_i^{\prime}+f_{1_i}^{\prime}+\lambda_1 \int_a^b\left(k_{1 i}^{\prime} u(t)+k_{2 i}^{\prime} u(t) v(t)\right) d t . \\ \boldsymbol{v}_{\boldsymbol{i}}^{(\mathbf{3})}=-q_{2_i}^{\prime} u_i^{\prime}-q_{2_i} u_i^{\prime \prime}-p_{2_i}^{\prime} v_i-p_{2_i} v_i^{\prime}+f_{2_i}^{\prime}+\lambda_2 \int_a^b\left(k_{3 i}^{\prime} u(t) u(t)+k_{4 i}^{\prime} u(t) v(t)\right) d t .\end{gathered}$           (40)

$\begin{aligned}& \boldsymbol{u}_{\boldsymbol{i}}^{(\mathbf{4})}=\left(p_{2_i} q_{1_i}-q_{1_i}^{\prime \prime}\right) v_i^{\prime}-2 q_{1_i}^{\prime} v_i^{\prime \prime}+q_{1_i} p_{2_i}^{\prime} v_i+\left(-p_{1_i}+q_{1_i} q_{2_i}\right) u_i^{\prime \prime}+\left(q_{1_i} q_{2_i}^{\prime}-2 p_{1_i}^{\prime}\right) u_i^{\prime}-p_{1_i}^{\prime \prime} u_i \\& +f_1^{\prime \prime}{ }_i-q_{1 i} f_2^{\prime}{ }_i-\lambda_2 q_{1_i} \int_a^b\left(k_{3 i}^{\prime} u(t)+k_{4 i}^{\prime} v(t)\right) d t+\lambda_1 \int_a^b\left(k_{1 i}^{\prime \prime} u(t)+k_{2 i}^{\prime \prime} v(t)\right) d t \\& \boldsymbol{v}_{\boldsymbol{i}}^{(\mathbf{4})}=\left(p_{1_i} q_{2_i}-q_{2_i}^{\prime \prime}\right) u_i^{\prime}-2 q_{2_i}^{\prime} u_i^{\prime \prime}+q_{2_i} p_{1_i}^{\prime} u_i+\left(-p_{2_i}+q_{1_i} q_{2_i}\right) v_i^{\prime \prime}+\left(q_{2_i} q_{1_i}^{\prime}-2 p_{2_i}^{\prime}\right) v_i^{\prime}-p_{2_i}^{\prime \prime} v_i \\& +f_2^{\prime \prime}{ }_i-q_{2_i} f_1^{\prime}{ }_i-\lambda_1 q_{2_i} \int_a^b\left(k_{1 i}^{\prime} u(t)+k_{2 i}^{\prime} v(t)\right) d t+\lambda_2 \int_a^b\left(k_{1 i}^{\prime \prime} u(t)+k_{2 i}^{\prime \prime} v(t)\right) d t .\end{aligned}$            (41)

Then substituting $\boldsymbol{u}_i^{(\mathbf{3})}$ and $\boldsymbol{v}_i^{(3)}$ from (40) and $\boldsymbol{u}_i^{(4)}$ and $\boldsymbol{v}_i^{(4)}$ from (41) into the term E4 in (38) and simplifying produces

$\begin{gathered}E_{1_4}=-\frac{h^2\left(p_{2_i} q_{1_i}+q_{1_i}^{\prime \prime}\right)}{12} v_i^{\prime}-\frac{h^2 q_{1_i}^{\prime}}{6} v_i^{\prime \prime}-\frac{h^2 q_{1_i} p_{2_i}^{\prime}}{12} v_i-\frac{h^2\left(p_{1_i}+q_{1_i} q_{2_i}\right)}{12} u_i^{\prime \prime} \\ -\frac{h^2\left(q_{1_i} q_{2_i}^{\prime}+2 p_{1_i}^{\prime}\right)}{12} u_i^{\prime}-\frac{h^2 p_{1_i}^{\prime \prime}}{12} u_i +\frac{h^2}{12} f_1^{\prime \prime}+\frac{h^2 q_{1_i}}{12} f_{2_i}^{\prime} \\ +\frac{\lambda_2 h^2 q_{1_i}}{12} \int_a^b\left(k_{3 i}^{\prime} u(t)+k_{4 i}^{\prime} v(t)\right) d t+\frac{\lambda_1 h^2}{12} \int_a^b\left(k_{1 i}^{\prime \prime} u(t)+k_{2 i}^{\prime \prime} v(t)\right) d t\end{gathered}$

$\begin{aligned} E_{24} & =-\frac{h^2\left(p_{1_i} q_{2_i}+q_{2_i}^{\prime \prime}\right)}{12} u_i^{\prime}-\frac{h^2 q_{2_i}^{\prime}}{6} u_i^{\prime \prime}-\frac{h^2 q_{2_i} p_{1_i}^{\prime}}{12} u_i \\ & -\frac{h^2\left(p_{2_i}+q_{1_i} q_{2_i}\right)}{12} v_i^{\prime \prime}-\frac{h^2\left(q_{2_i} q_{1 i}^{\prime}+2 p_{2_i}^{\prime}\right)}{12} v_i^{\prime} \\ & -\frac{h^2 p_{1_i}^{\prime \prime}}{12} v_i+\frac{h^2}{12} f_2^{\prime \prime}+\frac{h^2 q_{1_i}}{12} f_{1_i}^{\prime} \\ & +\frac{\lambda_1 h^2 q_{2_i}}{12} \int_a^b\left(k_{1 i}^{\prime} u(t)+k_{2 i}^{\prime} v(t)\right) d t \\ & +\frac{\lambda_2 h^2}{12} \int_a^b\left(k_{3 i}^{\prime \prime} u(t)+k_{4 i}^{\prime \prime} v(t)\right) d t\end{aligned}$

Finally, replace the terms $u_i^{\prime \prime}, u_i^{\prime}, v_i^{\prime \prime}$ and $v_i^{\prime}$ in $E_{1_4}$ and $E_{2_4}$ with those in (7) and (8) to obtain

$\begin{aligned} & \boldsymbol{E}_{1_4}=-\frac{h^2\left(p_{2_i} q_{1_i}+q_{1_i}^{\prime \prime}\right)}{12} \delta_x v_i-\frac{h^2 q_{1_i}^{\prime}}{6} \delta_x^2 v_i-\frac{h^2 q_{1_i} p_{2_i}^{\prime}}{12} v_i-\frac{h^2\left(p_{1_i}+q_{1_i} q_{2_i}\right)}{12} \delta_x^2 u_i-\frac{h^2\left(q_{1_i} q_{2_i}^{\prime}+2 p_{1_i}^{\prime}\right)}{12} \delta_x u_i \\ & -\frac{h^2 p_{1_i}^{\prime \prime}}{12} u_i+\frac{h^2}{12} f_1^{\prime \prime}+\frac{h^2 q_{1_i}}{12} f_{2_i}^{\prime}+\frac{\lambda_2 h^2 q_{1_i}}{12} \int_a^b\left(k_{3 i}^{\prime} u(t)+k_4^{\prime} v(t)\right) d t+\frac{\lambda_1 h^2}{12} \int_a^b\left(k_{1 i}^{\prime \prime} u(t)+k_{2 i}^{\prime \prime} v(t)\right) d t-\varepsilon r r 3 \end{aligned}$

$\begin{gathered}\boldsymbol{E}_{2_4}=-\frac{h^2\left(p_{1_i} q_{2_i}+q_{2_i}^{\prime \prime}\right)}{12} \delta_x u_i-\frac{h^2 q_{2_i}^{\prime}}{6} \delta_x^2 u_i-\frac{h^2 q_{2_i} p_{1_i}^{\prime}}{12} u_i-\frac{h^2\left(p_{2_i}+q_{1_i} q_{2_i}\right)}{12} \delta_x^2 v_i-\frac{h^2\left(q_{2_i} q_{1_i}^{\prime}+2 p_{2_i}^{\prime}\right)}{12} \delta_x v_i \\ -\frac{h^2 p_{1_i}^{\prime \prime}}{12} v_i+\frac{h^2}{12} f_{2_i^{\prime \prime}}^{\prime \prime}+\frac{h^2 q_{1_i}}{12} f_{1_i}^{\prime}+\frac{\lambda_1 h^2 q_{2_i}}{12} \int_a^b \int_a^b\left(k_{1 i}^{\prime} u(t)+k_{2 i}^{\prime} v(t)\right) d t+\frac{\lambda_2 h^2}{12} \int_a^b\left(k_{3 i}^{\prime \prime} u(t)+k_{4 i}^{\prime \prime} v(t)\right) d t-\varepsilon r r 4 \end{gathered}$

the error term Err3 and Err4 are given as

$\begin{gathered}\varepsilon \boldsymbol{r r} 3=\frac{h^2\left(p_{2_i} q_{1_i}+q_{1_i}^{\prime \prime}\right)}{12}\left(\frac{h^2}{3!} v_i^{(3)}+\frac{h^4}{5!} v_i^{(5)}+\frac{h^6}{7!} v_i^{(7)}+\frac{h^8}{9!} v_i^{(9)}+\frac{h^{10}}{11!} v_i^{(11)}+\cdots\right)+\frac{4 h^2 q_{1_i}^{\prime}}{6}\left(\frac{h^2}{4!} v_i^{(4)}+\frac{h^4}{6!} v_i^{(6)}+\frac{h^6}{8!} v_i^{(8)}+\frac{h^8}{10!} v_i^{(10)}+\frac{h^{10}}{12!} v_i^{(12)} \cdots\right) \\ +\frac{4 h^2\left(p_{1_i}+q_{1_i} q_{2_i}\right)}{12}\left(\frac{h^2}{4!} u_i^{(4)}+\frac{h^4}{6!} u_i^{(6)}+\frac{h^6}{8!} u_i^{(8)}+\frac{h^8}{10!} u_i^{(10)}+\frac{h^{10}}{12!} u_i^{(12)} \cdots\right)+\frac{h^2\left(q_{1_i} q_{2_i}^{\prime}+2 p_{1_i}^{\prime}\right)}{12}\left(\frac{h^2}{3!} u_i^{(3)}+\frac{h^4}{5!} u_i^{(5)}+\frac{h^6}{7!} u_i^{(7)}+\frac{h^8}{9!} u_i^{(9)}+\frac{h^{10}}{11!} u_i^{(11)}+\cdots\right)\end{gathered}$

$\begin{gathered}\varepsilon \boldsymbol{r r} \mathbf{4}=\frac{h^2\left(p_{1_i} q_{2_i}+q_{2_i}^{\prime \prime}\right)}{12}\left(\frac{h^2}{3!} u_i^{(3)}+\frac{h^4}{5!} u_i^{(5)}+\frac{h^6}{7!} u_i^{(7)}+\frac{h^8}{9!} u_i^{(9)}+\frac{h^{10}}{11!} u_i^{(11)}+\cdots\right)+\frac{2 h^2 q_{2_i}^{\prime}}{6}\left(\frac{h^2}{4!} u_i^{(4)}+\frac{h^4}{6!} u_i^{(6)}+\frac{h^6}{8!} u_i^{(8)}+\frac{h^8}{10!} u_i^{(10)}+\frac{h^{10}}{12!} u_i^{(12)} \ldots\right) \\ +\frac{2 h^2\left(p_{2_i}+q_{1_i} q_{2_i}\right)}{12}\left(\frac{h^2}{4!} v_i^{(4)}+\frac{h^4}{6!} v_i^{(6)}+\frac{h^6}{8!} v_i^{(8)}+\frac{h^8}{10!} v_i^{(10)}+\frac{h^{10}}{12!} v_i^{(12)} \cdots\right)+\frac{h^2\left(q_{2_i} q_{1_i}^{\prime}+2 p_{2_i}^{\prime}\right)}{12}\left(\frac{h^2}{3!} v_i^{(3)}+\frac{h^4}{5!} v_i^{(5)}+\frac{h^6}{7!} v_i^{(7)}+\frac{h^8}{9!} v_i^{(9)}+\frac{h^{10}}{11!} v_i^{(11)}+\cdots\right) .\end{gathered}$

Adding $\boldsymbol{E}_{\mathbf{1}_{\mathbf{4}}}$ and $\boldsymbol{E}_{\mathbf{2}_{\mathbf{4}}}$ to the second order central difference scheme, gives

$\begin{aligned} & \left(q_{1_i}+\frac{h^2\left(p_{2_i} q_{1_i}+q_{1_i}^{\prime \prime}\right)}{12}\right) \delta_x v_i-\frac{h^2 q_{1_i}^{\prime}}{6} \delta_x^2 v_i-\frac{h^2 q_{1_i} p_{2_i}^{\prime}}{12} v_i+\left(1+\frac{h^2\left(p_{1_i}+q_{1_i} q_{2_i}\right)}{12}\right) \delta_x^2 u_i+\frac{h^2\left(q_{1_i} q_{2_i}^{\prime}+2 p_{1_i}^{\prime}\right)}{12} \delta_x u_i \\ & +\left(1+\frac{h^2 p_{1_i}^{\prime \prime}}{12}\right) u_i+\frac{h^2}{12} f_1^{\prime \prime}+\frac{h^2 q_{1_i}}{12} f_{2_i}^{\prime}+\frac{\lambda_2 h^2 q_{1_i}}{12} \int_a^b\left(k_{3 i}^{\prime} u(t)+k_{4 i}^{\prime} v(t)\right) d t+\frac{\lambda_1 h^2}{12} \int_a^b\left(k_{1 i}^{\prime \prime} u(t)+k_{2 i}^{\prime \prime} v(t)\right) d t-\varepsilon r r 5\end{aligned}$

$\begin{aligned} &=f_{1_i}+\lambda_1 \int_a^b\left(k_{1_i}(t) u(t)+k_{2_i}(t) v(t)\right) d t,\left(q_{2_i}+\frac{h^2\left(p_{1_i} q_{2_i}+q_{2_i}^{\prime \prime}\right)}{12}\right) \delta_x u_i-\frac{h^2 q_{2_i}^{\prime}}{6} \delta_x^2 u_i-\frac{h^2 q_{2_i} p_{1 i}^{\prime}}{12} u_i+\left(1+\frac{h^2\left(p_{2_i}+q_{1_i} q_{2_i}\right)}{12}\right) \delta_x^2 v_i+\frac{h^2\left(q_{2_i} q_{1_i}^{\prime}+2 p_{2_i}^{\prime}\right)}{12} \delta_x v_i \\ &+\left(p_{2_i}+\frac{h^2 p_{1_i}^{\prime \prime}}{12}\right) v_i+\frac{h^2}{12} f_2^{\prime \prime}{ }_i+\frac{h^2 q_{1 i}}{12} f_{1_i}^{\prime}+\frac{\lambda_1 h^2 q_{2_i}}{12} \int_a^b\left(k_{1 i}^{\prime} u(t)+k_{2 i}^{\prime} v(t)\right) d t+\frac{\lambda_2 h^2}{12} \int_a^b\left(k_{3 i}^{\prime \prime} u(t)+k_{4 i}^{\prime \prime} v(t)\right) d t-\varepsilon r r 6=f_{2_i}+\lambda_2 \int_a^b\left(k_{3 i}(t) u(t)+k_{4 i}(t) v(t)\right) d t \end{aligned}$     (42)

If we replace the term $\delta_x u_i, \delta_x^2 u_i, \delta_x v_i$ and $\delta_x^2 v_i$ in the Eq. (36) by the terms in (4), and neglecting the error terms ($\varepsilon r r 1$) and ($\varepsilon r r 2$), the subsequent accurate central difference of second order schemes are obtained

$\begin{gathered}A_i u_{i+1}+B_i u_i+C_i u_{i-1}+D_i v_{i+1}+E_i v_{i-1}=F_{2_i}+I_{1_i} \\ G_i v_{i+1}+H_i v_i+I_i v_{i-1}+J_i u_{i+1}+K_i u_{i-1}=F_{2 i}+I_{2 i} \\ \mathrm{i}=1,2,3, \ldots, N\end{gathered}$    (43)

where the coefficients $A_i, B_i, C_i, F_i$ and $I_i$ and the following fourth-order compact finite difference scheme is defined as

$\begin{gathered}A_i=\frac{1}{h^2}+\frac{\left(p_{1_i}+q_{1_i} q_{2_i}\right)}{12}+\frac{h\left(q_{1_i} q_{2_i}^{\prime}+2 p_{1_i}^{\prime}\right)}{24}, B_i=\frac{-2}{h^2}-\frac{\left(p_{1_i}+q_{1_i} q_{2_i}\right)}{6}+p_{1_i}+\frac{h^2 p_{1_i}^{\prime \prime}}{12} \\ C_i=\frac{1}{h^2}+\frac{\left(p_{1_i}+q_{1_i} q_{2_i}\right)}{12}-\frac{h\left(q_{1_i} q_{2_i}^{\prime}+2 p_{1_i}^{\prime}\right)}{24}, D_i=-\frac{q_{1_i}^{\prime}}{6}+\frac{q_{1_i}}{2 h}+\frac{h\left(p_{2_i} q_{1_i}+q_{1_i}^{\prime \prime}\right)}{24}, E_i=\frac{q_{1_i}^{\prime}}{3}-\frac{h^2 q_{1_i} p_{2_i}^{\prime}}{12} \\ L_i=-\frac{q_{1_i}^{\prime}}{6}-\frac{q_{1_i}}{2 h}-\frac{h\left(p_{2_i} q_{1_i}+q_{1_i}^{\prime \prime}\right)}{24}, F_{1_i}=f_{1_i}-\frac{h^2}{12} f_{1_i}^{\prime \prime}-\frac{h^2 q_{1_i}}{12} f_{2_i}^{\prime}\end{gathered}$

$\begin{gathered}I_{1_i}=\lambda_1 \int_a^b\left(k_{1_i}(t) u(t)+k_{2_i}(t) v(t)\right) d t-\frac{\lambda_2 h^2 q_{1_i}}{12} \int_a^b\left(k_{3 i}^{\prime} u(t)+k_{4 i}^{\prime} v(t)\right) d t-\frac{\lambda_1 h^2}{12} \int_a^b\left(k_{1 i}^{\prime \prime} u(t)+k_{2 i}^{\prime \prime} v(t)\right) d t \\ G_i=\frac{1}{h^2}+\frac{\left(p_{2_i}+q_{1_i} q_{2_i}\right)}{12}+\frac{h\left(q_{2_i} q_{1_i}^{\prime}+2 p_{2_i}^{\prime}\right)}{24}, H_i=\frac{-2}{h^2}-\frac{\left(p_{2_i}+q_{1_i} q_{2_i}\right)}{6}+p_{2_i}+\frac{h^2 p_{1_i}^{\prime \prime}}{12} \\ I_i=\frac{1}{h^2}+\frac{\left(p_{2_i}+q_{1_i} q_{2_i}\right)}{12}-\frac{h\left(q_{2_i} q_{1_i}^{\prime}+2 p_{2_i}^{\prime}\right)}{24}, J_i=-\frac{q_{2_i}^{\prime}}{6}+\frac{q_{2_i}}{2 h}+\frac{h\left(p_{1_i} q_{2_i}+q_{2_i}^{\prime \prime}\right)}{24}\end{gathered}$

$\begin{gathered}K_i=\frac{h^2 q_{2_i} p_{1_i}^{\prime}}{12}+\frac{q_{2_i}^{\prime}}{3}, M_i=-\frac{q_{2_i}^{\prime}}{6}-\frac{q_{2_i}}{2 h}-\frac{h\left(p_{1_i} q_{2_i}+q_{2_i}^{\prime \prime}\right)}{24}, F_{2_i}=f_{2_i}-\frac{h^2}{12} f_{2_i}^{\prime \prime}-\frac{h^2 q_{1_i}}{12} f_{1_i}^{\prime} \\ I_{2_i}=\lambda_2 \int_a^b\left(k_{3_i}(t) u(t)+k_{4_i}(t) v(t)\right) d t+\frac{\lambda_1 h^2 q_{2_i}}{12} \int_a^b\left(k_{1 i}^{\prime} u(t)+k_{2 i}^{\prime} v(t)\right) d t+\frac{\lambda_2 h^2}{12} \int_a^b\left(k_{3 i}^{\prime \prime} u(t)+k_{4 i}^{\prime \prime} v(t)\right) d t\end{gathered}$     (44)

and the new error term Err 5 and Err 6 is a linear combination of the error obtained from the central difference scheme (11) and the error obtained from E4 i.e. $\operatorname{Err} 5=\operatorname{Err} 1+\operatorname{Err} 3$ and $\operatorname{Err} 6=\operatorname{Err} 2+\operatorname{Err} 4$ so that

$\begin{aligned} \operatorname{Err} \mathbf{5}=h^4 & \left(\frac{2\left(3 q_{1 i} v_i^{(5)}+u_i^{(6)}\right)}{6!}+\frac{k_1 u_i^{(3)}+k_2 v_i^{(3)}}{3!}\right. \\ & \left.+\frac{2\left(k_3 u_i^{(4)}+k_4 v_i^{(4)}\right)}{4!}\right) \\ & +h^6\left(\frac{2\left(4 q_{1 i} v_i^{(7)}+u_i^{(8)}\right)}{8!}\right. \\ & \left.+\frac{k_1 u_i^{(5)}+k_2 v_i^{(5)}}{5!}+\frac{2\left(k_3 u_i^{(6)}+k_4 v_i^{(6)}\right)}{6!}\right) \\ & +h^8\left(\frac{2\left(5 q_1 v_i^{(9)}+u_i^{(10)}\right)}{10!}\right. \\ & \left.+\frac{k_1 u_i^{(7)}+k_2 v_i^{(7)}}{7!}+\frac{2\left(k_3 u_i^{(8)}+k_4 v_i^{(8)}\right)}{8!}\right) \\ & +h^{10}\left(\frac{2\left(6 q_{1 i} v_i^{(11)}+u_i^{(12)}\right)}{12!}\right. \\ & +\frac{k_1 u_i^{(9)}+k_2 v_i^{(9)}}{10!} \\ & \left.+\frac{2\left(k_3 u_i^{(10)}+k_4 v_i^{(10)}\right)}{10!}\right)+\cdots\end{aligned}$

$\begin{aligned} & \boldsymbol{\operatorname { E r r }} \mathbf{6}=h^4\left(\frac{2\left(3 q_{2_i} u_i^{(5)}+v_i^{(6)}\right)}{6!}+\frac{k_5 u_i^{(3)}+k_6 v_i^{(3)}}{3!}+\frac{2\left(k_7 u_i^{(4)}+k_8 v_i^{(4)}\right)}{4!}\right) \\ & +h^6\left(\frac{2\left(4 q_{1_i} u_i^{(7)}+v_i^{(8)}\right)}{8!}+\frac{k_5 u_i^{(5)}+k_6 v_i^{(5)}}{5!}+\frac{2\left(k_7 u_i^{(6)}+k_8 v_i^{(6)}\right)}{6!}\right) \\ & + \\ & +h^8\left(\frac{2\left(5 q_{1 i} u_i^{(9)}+v_i^{(10)}\right)}{10!}+\frac{k_5 u_i^{(7)}+k_6 v_i^{(7)}}{7!}+\frac{2\left(k_7 u_i^{(8)}+k_8 v_i^{(8)}\right)}{8!}\right) \\ & +h^{10}\left(\frac{2\left(6 q_{1_i} u_i^{(11)}+u_i^{(12)}\right)}{12!}+\frac{k_5 u_i^{(9)}+k_6 v_i^{(9)}}{10!}+\frac{2\left(k_7 u_i^{(10)}+k_8 v_i^{(10)}\right)}{10!}\right)+\cdots\end{aligned}$

where $k_1=\frac{q_{1 i} q_{2_i}^{\prime}+2 p_{1_i}^{\prime}}{12}, k_2=\frac{p_{2_i} q_{1_i}+q_{1 i}^{\prime \prime}}{12}, k_3=\frac{2\left(p_{1_i}+q_{1_i} q_{2_i}\right)}{12}, k_4=\frac{2 q_{1 i}^{\prime}}{6}, k_5=\frac{p_{1_i} q_{2_i}+q_{2 i}^{\prime \prime}}{12}, k_6=\frac{q_{2_i} q_{1_i^{\prime}}^{\prime}+2 p_{2_i}^{\prime}}{12}, k_7=\frac{2 q_{2_i}^{\prime}}{6}$, and $k_8=$ $\frac{2\left(p_{2_i}+q_{1_i} q_{2_i}\right)}{12}$,

Here, $u_0, v_0, u_n$ and $v_n$ are Dirichlet boundary conditions in general system, when $i=1$, and $i=n-1$ we obtain $u_0$ and $u_n$ respectively. The system in Eqs. (35) and (36) consists of ($2N-2$) linear equations with ($2N-2$) unknowns $\left(u_1, u_2, u_3, \ldots, u_{n-2}, u_{n-1}, v_1, v_2, v_3, \ldots, v_{n-2}, v_{n-1}\right)$ and can be written in the following matrix form

$\left[\begin{array}{ll}\gamma^1+\mathcal{A} & \delta^1+\mathcal{B} \\ \gamma^2+\mathcal{C} & \delta^2+\mathcal{D}\end{array}\right]\left[\begin{array}{l}U \\ V\end{array}\right]=\left[\begin{array}{l}\mathcal{F}_1+\omega_1+\varphi_1 \\ \mathcal{F}_2+\omega_2+\varphi_2\end{array}\right]$,

where,

$\gamma^1=\left[\begin{array}{ccccc}\gamma_2 & \gamma_3 & 0 & \cdots & 0 \\ \gamma_1 & \gamma_2 & \gamma_3 & \cdots & 0 \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & \cdots & \gamma_1 & \gamma_2 & \gamma_3 \\ 0 & \cdots & 0 & \gamma_1 & \gamma_2\end{array}\right], \delta^1=\left[\begin{array}{ccccc}\delta_2 & \delta_3 & 0 & \cdots & 0 \\ \delta_1 & \delta_2 & \delta_3 & \cdots & 0 \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & \cdots & \delta_1 & \delta_2 & \delta_3 \\ 0 & \cdots & 0 & \delta_1 & \delta_2\end{array}\right]$,

$\mathcal{A}=\left[\begin{array}{ccccccc}a_{1,1} & b_{1,2} & a_{1,3} & b_{1,4} & \cdots & b_{1, n-2} & a_{1, n-1} \\ a_{2,1} & b_{1,2} & a_{3,3} & b_{2,4} & \cdots & b_{2, n-2} & a_{2, n-1} \\ a_{3,1} & b_{1,2} & a_{3,3} & b_{3,4} & \cdots & b_{3, n-2} & a_{3, n-1} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ a_{n-3,1} & b_{n-3,2} & a_{n-3,3} & b_{n-3,4} & \cdots & b_{n-3, n-2} & a_{n-3, n-1} \\ a_{n-2,1} & b_{n-2,2} & a_{n-2,3} & b_{n-2,4} & \cdots & b_{n-2, n-2} & a_{n-2, n-1} \\ a_{n-1,1} & b_{n-1,2} & a_{n-1,3} & b_{n-1,4} & \cdots & b_{n-1, n-2} & a_{n-1, n-1}\end{array}\right]$,

$\mathcal{B}=\left[\begin{array}{ccccccc}c_{1,1} & d_{1,2} & c_{1,3} & d_{1,4} & \cdots & d_{1, n-2} & c_{1, n-1} \\ c_{2,1} & d_{1,2} & c_{3,3} & d_{2,4} & \cdots & d_{2, n-2} & c_{2, n-1} \\ c_{3,1} & d_{1,2} & c_{3,3} & d_{3,4} & \cdots & d_{3, n-2} & c_{3, n-1} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ c_{n-3,1} & d_{n-3,2} & c_{n-3,3} & d_{n-3,4} & \cdots & d_{n-3, n-2} & c_{n-3, n-1} \\ c_{n-2,1} & d_{n-2,2} & c_{n-2,3} & d_{n-2,4} & \cdots & d_{n-2, n-2} & c_{n-2, n-1} \\ c_{n-1,1} & d_{n-1,2} & c_{n-1,3} & d_{n-1,4} & \cdots & d_{n-1, n-2} & c_{n-1, n-1}\end{array}\right]$,

$U=\left[\begin{array}{c}u_1 \\ u_2 \\ u_3 \\ \vdots \\ u_{n-3} \\ u_{n-2} \\ u_{n-1}\end{array}\right], V=\left[\begin{array}{c}v_1 \\ v_2 \\ v_3 \\ \vdots \\ v_{n-3} \\ v_{n-2} \\ v_{n-1}\end{array}\right], \mathcal{F}_1=\left[\begin{array}{c}f_{1_1} \\ f_{1_2} \\ f_{1_3} \\ \vdots \\ f_{1_{n-3}} \\ f_{1_{n-2}} \\ f_{1_{n-1}}\end{array}\right], \mathcal{F}_2=\left[\begin{array}{c}f_{2_1} \\ f_{2_2} \\ f_{2_3} \\ \vdots \\ f_{2_{n-3}} \\ f_{2_{n-2}} \\ f_{2_{n-1}}\end{array}\right]$,

$\omega_1=\left[\begin{array}{c}\left(\alpha_{1_1}-\gamma_1\right) u_0+\beta_{1_1} u_n \\ \alpha_{1_2} u_0+\beta_{1_2} u_n \\ \alpha_{1_3} u_0+\beta_{1_3} u_n \\ \vdots \\ \alpha_{1_{n-3}} u_0+\beta_{1_{n-3}} u_n \\ \alpha_{1_{n-2}} u_0+\beta_{1_{n-2}} u_n \\ \alpha_{1_{n-1}} u_0+\left(\beta_{1_{n-1}}-\gamma_3\right) u_n\end{array}\right], \varphi_1=\left[\begin{array}{c}\left(\alpha_{2_1}-\delta_1\right) v_0+\beta_{2_1} v_n \\ \alpha_{2_2} v_0+\beta_{2_2} v_n \\ \alpha_{2_3} v_0+\beta_{2_3} v_n \\ \vdots \\ \alpha_{2_{n-3}} v_0+\beta_{2_{n-3}} v_n \\ \alpha_{2_{n-2}} v_0+\beta_{2_{n-2}} v_n \\ \alpha_{2_{n-1}} v_0+\left(\beta_{2_{n-1}}-\delta_3\right) v_n\end{array}\right]$,

$\gamma^2=\left[\begin{array}{ccccc}\gamma_5 & \gamma_6 & 0 & \cdots & 0 \\ \gamma_4 & \gamma_5 & \gamma_6 & \cdots & 0 \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & \cdots & \gamma_4 & \gamma_5 & \gamma_6 \\ 0 & \cdots & 0 & \gamma_4 & \gamma_5\end{array}\right], \omega_2=\left[\begin{array}{c}\left(\alpha_{3_1}-\gamma_4\right) u_0+\beta_{3_1} u_n \\ \alpha_{3_2} u_0+\beta_{3_2} u_n \\ \alpha_{3_3} u_0+\beta_{3_3} u_n \\ \vdots \\ \alpha_{3_{n-3}} u_0+\beta_{3_{n-3}} u_n \\ \alpha_{3_{n-2}} u_0+\beta_{3_{n-2}} u_n \\ \alpha_{3_{n-1}} u_0+\left(\beta_{3_{n-1}}-\gamma_6\right) u_n\end{array}\right]$,

$\mathcal{C}=\left[\begin{array}{ccccccc}r_{1,1} & s_{1,2} & r_{1,3} & s_{1,4} & \cdots & s_{1, n-2} & r_{1, n-1} \\ r_{2,1} & s_{1,2} & r_{3,3} & s_{2,4} & \cdots & s_{2, n-2} & r_{2, n-1} \\ r_{3,1} & s_{1,2} & r_{3,3} & s_{3,4} & \cdots & s_{3, n-2} & r_{3, n-1} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ r_{n-3,1} & s_{n-3,2} & r_{n-3,3} & s_{n-3,4} & \cdots & s_{n-3, n-2} & r_{n-3, n-1} \\ r_{n-2,1} & s_{n-2,2} & r_{n-2,3} & s_{n-2,4} & \cdots & s_{n-2, n-2} & r_{n-2, n-1} \\ r_{n-1,1} & s_{n-1,2} & r_{n-1,3} & s_{n-1,4} & \cdots & s_{n-1, n-2} & r_{n-1, n-1}\end{array}\right]$,

$\varphi_2=\left[\begin{array}{c}\left(\alpha_{41}-\delta_4\right) v_0+\beta_{41} v_n \\ \alpha_{42} v_0+\beta_{42} v_n \\ \alpha_{43} v_0+\beta_{43} v_n \\ \vdots \\ \alpha_{4 n-3} v_0+\beta_{4 n-3} v_n \\ \alpha_{4 n-2} v_0+\beta_{4 n-2} v_n \\ \alpha_{4 n-1} v_0+\left(\beta_{4 n-1}-\delta_6\right) v_n\end{array}\right], \delta^2=\left[\begin{array}{ccccc}\delta_5 & \delta_6 & 0 & \cdots & 0 \\ \delta_4 & \delta_5 & \delta_6 & \cdots & 0 \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & \cdots & \delta_4 & \delta_5 & \delta_6 \\ 0 & \cdots & 0 & \delta_4 & \delta_5\end{array}\right]$,

$\mathcal{D}=\left[\begin{array}{ccccccc}y_{1,1} & z_{1,2} & y_{1,3} & z_{1,4} & \cdots & z_{1, n-2} & y_{1, n-1} \\ y_{2,1} & z_{1,2} & y_{3,3} & z_{2,4} & \cdots & z_{2, n-2} & y_{2, n-1} \\ y_{3,1} & z_{1,2} & y_{3,3} & z_{3,4} & \cdots & z_{3, n-2} & y_{3, n-1} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ y_{n-3,1} & z_{n-3,2} & y_{n-3,3} & z_{n-3,4} & \cdots & z_{n-3, n-2} & y_{n-3, n-1} \\ y_{n-2,1} & z_{n-2,2} & y_{n-2,3} & z_{n-2,4} & \cdots & z_{n-2, n-2} & y_{n-2, n-1} \\ y_{n-1,1} & z_{n-1,2} & y_{n-1,3} & z_{n-1,4} & \cdots & z_{n-1, n-2} & y_{n-1, n-1}\end{array}\right]$,

$\begin{gathered}a_{i, j}=\frac{h^3 \lambda_1}{9}\left(k_{1_{i+1,2 j-1}}+k_{1_{i-1,2 j-1}}\right)+\frac{10 h^3 \lambda_1}{9} k_{1_{i, 2 j-1}}+\frac{q_{1_i} h^4 \lambda_2}{18}\left(k_{3_{i+1,2 j-1}}+k_{3_{i-1,2 j-1}}\right), \\ b_{i, j}=\frac{h^3 \lambda_1}{18}\left(k_{1_{i+1,2 j}}+k_{1_{i-1,2 j}}\right)+\frac{10 h^3 \lambda_1}{18} k_{1_{i, 2 j}}+\frac{q_{1_i} h^4 \lambda_2}{36}\left(k_{3_{i+1,2 j}}+k_{3_{i-1,2 j}}\right),\end{gathered}$

$\begin{gathered}c_{i, j}=\frac{h^3 \lambda_1}{9}\left(k_{2_{i+1,2 j-1}}+k_{2_{i-1,2 j-1}}\right)+\frac{10 h^3 \lambda_1}{9} k_{2_{i, 2 j-1}}+\frac{q_{1_i} h^4 \lambda_2}{18}\left(k_{4_{i+1,2 j-1}}+k_{4_{i-1,2 j-1}}\right), \\ d_{i, j}=\frac{h^3 \lambda_1}{18}\left(k_{2_{i+1,2 j}}+k_{2_{i-1,2 j}}\right)+\frac{10 h^3 \lambda_1}{18} k_{2_{i, 2 j}}+\frac{q_{1_i} h^4 \lambda_2}{36}\left(k_{4_{i+1,2 j}}+k_{4_{i-1,2 j}}\right),\end{gathered}$

$\begin{gathered}w_{i, j}=\frac{h^3 \lambda_2}{9}\left(k_{3_{i+1,2 j-1}}+k_{3_{i-1,2 j-1}}\right)+\frac{10 h^3 \lambda_2}{9} k_{3_{i, 2 j-1}}+\frac{q_{2_i} h^4 \lambda_1}{18}\left(k_{1_{i+1,2 j-1}}+k_{1_{i-1,2 j-1}}\right), \\ x_{i, j}=\frac{h^3 \lambda_2}{18}\left(k_{3_{i+1,2 j}}+k_{3_{i-1,2 j}}\right)+\frac{10 h^3 \lambda_1}{18} k_{3_{i, 2 j}}+\frac{q_2 h^4 \lambda_1}{36}\left(k_{1_{i+1,2 j}}+k_{1_{i-1,2 j}}\right),\end{gathered}$

$\begin{gathered}y_{i, j}=\frac{h^3 \lambda_2}{9}\left(k_{4_{i+1,2 j-1}}+k_{4_{i-1,2 j-1}}\right)+\frac{10 h^3 \lambda_2}{9} k_{4_{i, 2 j-1}}+\frac{q_{2_i} h^4 \lambda_1}{18}\left(k_{2_{i+1,2 j-1}}+k_{2_{i-1,2 j-1}}\right), \\ z_{i, j}=\frac{h^3 \lambda_2}{18}\left(k_{4_{i+1,2 j}}+k_{4_{i-1,2 j}}\right)+\frac{10 h^3 \lambda_2}{18} k_{4_{i, 2 j}}+\frac{q_{2_i} h^4 \lambda_1}{36}\left(k_{2_{i+1,2 j}}+k_{2_{i-1,2 j}}\right),\end{gathered}$

$f_{1_i}=\frac{h^2}{12}\left(f_{1_{i+1}}+10 f_{1_i}+f_{1_{i-1}}\right)+\frac{q_{1_i} h^3}{24}\left(f_{2_{i+1}}-f_{2_{i-1}}\right), f_{2_i}=\frac{h^2}{12}\left(f_{2_{i+1}}+10 f_{2_i}+f_{2_{i-1}}\right)+\frac{q_{2_i} h^3}{24}\left(f_{1_{i+1}}-f_{1_{i-1}}\right)$,

$\begin{gathered}\boldsymbol{\gamma}_1=1+\frac{h^2\left(q_{1_i} q_{2_i}+p_{1_i}\right)}{12}+\frac{h^3\left(q_{1_i} q_{2_i}^{\prime}+2 p_{1_i}^{\prime}\right)}{24}, \gamma_2=-2-\frac{h^2\left(q_{1_i} q_{2_i}+p_{1_i}\right)}{6}+p_{1_i} h^2+\frac{h^4 p_{1_i}^{\prime \prime}}{12}, \\ \gamma_3=1+\frac{h^2\left(q_{1_i} q_{2_i}+p_{1_i}\right)}{12}-\frac{h^3\left(q_{1_i} q_{2_i}^{\prime}+2 p_{1_i}^{\prime}\right)}{24},\end{gathered}$

$\boldsymbol{\delta}_{\mathbf{1}}=\frac{h^2 q_{1_i}^{\prime}}{6}+\frac{h q_{1_i}}{2}+\frac{h^3\left(p_{2_i} q_{1_i}+q_{1_i}^{\prime \prime}\right)}{24}, \boldsymbol{\delta}_{\mathbf{2}}=\frac{-h^2 q_{1_i}^{\prime}}{3}+\frac{h^4 q_{1_i} p_{2_i}^{\prime}}{12}, \boldsymbol{\delta}_{\mathbf{3}}=\frac{h^2 q_{1_i}^{\prime}}{6}-\frac{h q_{1_i}}{2}-\frac{h^3\left(p_{2_i} q_{1_i}+q_{1_i}^{\prime \prime}\right)}{24}$,

$\boldsymbol{\gamma}_4=\frac{h^2 q_{2_i}^{\prime}}{6}+\frac{h q_{2_i}}{2}+\frac{h^3\left(p_{1_i} q_{2_i}+q_{2_i}^{\prime \prime}\right)}{24}, \boldsymbol{\gamma}_5=\frac{-h^2 q_{2_i}^{\prime}}{3}+\frac{h^4 q_{2_i} p_{1_i}^{\prime}}{12}, \boldsymbol{\gamma}_6=\frac{h^2 q_{2_i}^{\prime}}{6}-\frac{h q_{2_i}}{2}-\frac{h^3\left(p_{1_i} q_{2_i}+q_{2_i}^{\prime \prime}\right)}{24}$,

$\begin{gathered}\boldsymbol{\delta}_{\mathbf{4}}=1+\frac{h^2\left(q_{1_i} q_{2_i}+p_{2_i}\right)}{12}-\frac{h^3\left(q_{1_i}^{\prime} q_{2_i}+2 p_{2_i}^{\prime}\right)}{24}, \boldsymbol{\delta}_{\mathbf{5}}=-2-\frac{h^2\left(q_{1_i} q_{2_i}+p_{2_i}\right)}{6}+p_{2_i} h^2+\frac{h^4 p_{2_i}^{\prime \prime}}{12}, \\ \boldsymbol{\delta}_{\mathbf{6}}=1+\frac{h^2\left(q_{1_i} q_{2_i}+p_{2_i}\right)}{12}-\frac{h^3\left(q_{1_i}^{\prime} q_{2_i}+2 p_{2_i}^{\prime}\right)}{24}, \end{gathered}$

$\begin{gathered}\alpha_{1_i}=\frac{h^3 \lambda_1}{36}\left(k_{1_{i+1,0}}+k_{1_{i-1,0}}\right)+\frac{10 h^3 \lambda_1}{36} k_{1_{i, 0}} +\frac{q_{1 i} h^4 \lambda_2}{72}\left(k_{3_{i+1,0}}-k_{3_{i-1,0}}\right),\end{gathered}$

$\begin{gathered}\alpha_{2_i}=\frac{h^3 \lambda_{11}}{36}\left(k_{2_{i+1,0}}+k_{2_{i-1,0}}\right)+\frac{10 h^3 \lambda_1}{36} k_{2_{i, 0}} +\frac{q_{1 i} h^4 \lambda_2}{72}\left(k_{4_{i+1,0}}-k_{4_{i-1,0}}\right),\end{gathered}$

$\begin{gathered}\boldsymbol{\beta}_{1_i}=\frac{h^3 \lambda_1}{36}\left(k_{1_{i+1, n}}+k_{1_{i-1, n}}\right)+\frac{10 h^3 \lambda_1}{36} k_{1_{i, n}} +\frac{q_{1_i} h^4 \lambda_2}{72}\left(k_{3_{i+1, n}}-k_{3_{i-1, n}}\right),\end{gathered}$

$\begin{gathered}\boldsymbol{\beta}_{2_i}=\frac{h^3 \lambda_1}{36}\left(k_{2_{i+1, n}}+k_{2_{i-1, n}}\right)+\frac{10 h^3 \lambda_1}{36} k_{2_{i, n}} +\frac{q_{1_i} h^4 \lambda_2}{72}\left(k_{4_{i+1, n}}-k_{4_{i-1, n}}\right),\end{gathered}$

$\begin{gathered}\boldsymbol{a}_{3_i}=\frac{h^3 \lambda_2}{36}\left(k_{3_{i+1,0}}+k_{3_{i-1,0}}\right)+\frac{10 h^3 \lambda_2}{36} k_{3_{i, 0}} +\frac{q_{2_i} h^4 \lambda_1}{72}\left(k_{1_{i+1,0}}-k_{1_{i-1,0}}\right),\end{gathered}$

$\begin{gathered}\boldsymbol{a}_{4_i}=\frac{h^{3 \lambda_2}}{36}\left(k_{4_{i+1,0}}+k_{4_{i-1,0}}\right)+\frac{10 h^3 \lambda_2}{36} k_{4_{i, 0}} +\frac{q_{2_i} h^4 \lambda_1}{72}\left(k_{2_{i+1,0}}-k_{2_{i-1,0}}\right),\end{gathered}$

$\begin{gathered}\boldsymbol{\beta}_{3_i}=\frac{h^3 \lambda_2}{36}\left(k_{3_{i+1, n}}+k_{3_{i-1, n}}\right)+\frac{10 h^3 \lambda_2}{36} k_{3_{i, n}} +\frac{q_{2_i} h^4 \lambda_1}{72}\left(k_{1_{i+1, n}}-k_{1_{i-1, n}}\right),\end{gathered}$

$\begin{gathered}\boldsymbol{\beta}_{\mathbf{4}_{\boldsymbol{i}}}=\frac{h^3 \lambda_2}{36}\left(k_{4_{i+1, n}}+k_{4_{i-1, n}}\right)+\frac{10 h^{3 \lambda_2}}{36} k_{4_{i, n}}  +\frac{q_{2_i} h^4 \lambda_1}{72}\left(k_{2_{i+1, n}}-k_{2_{i-1, n}}\right) .\end{gathered}$

Algorithm 2

Input: $N, a, b$, and boundary conditions $u_0, u_n, v_0$ and $v_n$ where $\left(a=x_0, b=x_n\right)$.

Set: $h=\frac{b-a}{N}$.

for $i \longleftarrow 0$ to N do

     for $i \longleftarrow 0$ to N do     

$\begin{aligned} x_i & =a+i h \\ t_j & =a+j h\end{aligned}$

       end for

end for

for $i \longleftarrow 0$ to N do

$u_i \leftarrow u\left(x_i\right) ; v_i \leftarrow v\left(x_i\right)$

          for $i \longleftarrow 0$ to N do

$\begin{aligned} & k_{1_{i, j}} \leftarrow k_1\left(x_i, t_j\right) ; k_{2 i, j} \leftarrow k_2\left(x_i, t_j\right) ; \\ & k_{3_{i, j}} \leftarrow k_3\left(x_i, t_j\right) ; k_{4_{i, j}} \leftarrow k_4\left(x_i, t_j\right)\end{aligned}$

end for

end for

for $i \longleftarrow 0$ to N-1 do

     $\begin{aligned} & C_{1_i}=\mathcal{F}_1+\omega_1+\varphi_1 \\ & C_{2_i}=\mathcal{F}_2+\omega_2+\varphi_2\end{aligned}$

     for $i \longleftarrow 0$ to N-1 do

$\begin{gathered}\mathcal{A}=\gamma^1+\mathcal{A} \\ \mathcal{B}=\delta^1+\mathcal{B} \\ \mathcal{C}=\gamma^2+\mathcal{C} \\ \mathcal{D}=\delta^2+\mathcal{D}\end{gathered}$

     end for

end for

$\mathfrak{R}=[\mathcal{A}\ \mathcal{B} ; \mathcal{C}\ \mathcal{D}]$

$\begin{gathered}C_{1_i}=\left[C_{1_1}-\gamma_1 u_0-\delta_1 v_0 ; C_{1_2}: C_{1_{n-2}} ; C_{1_{n-1}}\right. \\ \left.-\gamma_3 u_n-\delta_3 v_n\right]\end{gathered}$

$\begin{gathered}C_{2_i}=\left[C_{2_1}-\gamma_4 u_0-\delta_4 v_0 ; C_{2_2}: C_{2_{n-2}} ; C_{2_{n-1}}\right. \\ \left.-\gamma_6 u_n-\delta_6 v_n\right]\end{gathered}$

$C=\left[C_{1_i} ; C_{2_i}\right]$

$\mathcal{W}=[U ; V]$

Output: $\mathcal{W} \leftarrow \Re \backslash C$

Next, the following examples, a system of Integro-Differential Equations, are solved by Compact Finite Difference Method.

Example 4: We deal with $p_1(x), p_2(x)=0, q_1(x)$ and $q_2(x) \neq 0$ such that:

$u^{\prime \prime}(x)+v^{\prime}(x)=2\left(e^x-\sin x\right)-\int_0^\pi e^x(u(t)-v(t)) d t$

$v^{\prime \prime}(x)+2 u^{\prime}(x)=\left(1+\frac{\pi}{2}\right) \cos x-\frac{\pi}{2} \sin x-\int_0^\pi \cos (x+t)(u(t)+v(t)) d t$,

$0 \leq x \leq \pi$

with Dirichlet boundary conditions,

$u(0)=0, \quad u(\pi)=0, \quad v(0)=1, \quad v(\pi)=-1$,

and the exact solutions are $u(x)=\sin x$ and $v(x)=\cos x$.

Example 5: We deal with $p_1(x), p_2(x)=0, q_1(x)$ and $q_2(x)=0$ such that:

$\begin{gathered}u^{\prime \prime}(x)=\frac{3 x}{10}+6-2 \int_0^1 x t(u(t)-3 v(t)) d t \\ v^{\prime \prime}(x)=15 x+\frac{4}{5}-3 \int_0^1\left(2 x+t^2\right)(u(t)-2 v(t)) d t \\ 0 \leq x \leq 1\end{gathered}$

with Dirichlet boundary conditions,

$u(0)=1, \quad u(1)=4, \quad v(0)=-1, \quad v(1)=2$,

and exact solutions are: $u(x)=3 x^2+1$ and $v(x)=x^3+2 x-1$.

Example 6: We deal with $p_1(x), p_2(x)=0, q_1(x)$ and $q_2(x) \neq 0$ such that:

$\begin{gathered}u^{\prime \prime}(x)+v^{\prime}(x)=e^x+\frac{x}{2}+\int_0^1 x t(u(t)+2 v(t)) d t, \\ v^{\prime \prime}(x)+u^{\prime}(x)=2-x e^1+\frac{x}{3}+e^x+\int_0^1 x(u(t)+2 v(t)) d t, \\ 0 \leq x \leq 1\end{gathered}$

with Dirichlet boundary conditions:

$u(0)=1, \quad u(1)=e^1, \quad v(0)=0, \quad v(1)=1$,

and the exact solutions are $u(x)=e^x$ and $v(x)=x^2$.

Table 1. Rate convergence for Example 1

$N$

$l^2(\mathrm{E}, h)$

Rate

$l^{\infty}(\mathrm{E}, \mathrm{h})$

Rate

10

2.5935e-07

 

3.5847e-07

 

20

1.6227e-08

3.9984

2.2488e-08

3.9946

40

1.0145e-09

3.9996

1.4079e-09

3.9975

80

6.3409e-11

3.9999

8.8014e-11

3.9997

Table 2. Rate convergence for Example 2

$N$

$l^2(\mathrm{E}, h)$

Rate

$l^{\infty}(\mathrm{E}, \mathrm{h})$

Rate

10

5.4065e-05

 

7.9096e-05

 

20

3.3780e-06

4.0005

5.0185e-06

3.9783

40

2.1112e-07

4.0000

3.1365e-07

4.0000

80

1.3195e-08

4

1.9603e-08

4.0000

Table 3. Rate convergence for Example 3

$N$

$l^2(\mathrm{E}, h)$

Rate

$l^{\infty}(\mathrm{E}, \mathrm{h})$

Rate

10

2.9359e-03

 

3.5700e-03

 

20

2.1150e-04

3.7951

2.4212e-04

3.8821

40

1.3614e-05

3.9575

1.6098e-05

3.9108

80

8.5649e-07

3.9905

1.0112e-06

3.9927

Table 4. Comparison of the pointwise error $\varepsilon(x)$ when $h=0.1$ and 0.05 for Example 1

$x_i$

Present Method

Reference [13]

Reference [18]

Reference [25]

N=10

N=20

N=10

N=20

N=10

N=10

2.05e-02

0.1

9.77e-08

6.11e-09

4.31e-05

1.21e-06

4.39e-06

0.2

1.86e-07

1.16e-08

3.07e-05

6.82e-07

1.81e-05

0.3

2.61e-07

1.63e-08

1.43e-05

1.41e-07

4.24e-05

0.4

3.18e-07

1.99e-08

2.37e-06

4.02e-07

7.83 e-05

0.5

3.52e-07

2.20e-08

1.91e-05

9.45e-07

1.27e-04

0.6

3.58e-07

2.24e-08

3.59e-05

1.49e-06

1.91e-04

0.7

3.31e-07

2.07e-08

5.28e-05

2.04e-06

2.70e-04

0.8

2.65e-07

1.66e-08

6.93e-05

2.59e-06

3.67 e-04

0.9

1.57e-07

9.82e-09

7.94e-05

3.12e-06

4.85e-04

Table 5. Comparison of the pointwise error $\varepsilon(x)$ when $h=0.0417$ and 0.0208 for Example 2

$x_i$

Present Method

Reference [13]

Reference [26]

Reference [26]

N=24

N=48

N=24

2.41e-06

 

N=48

1.51e-07

N=24

3.25e-06

 

N=48

2.41e-07

N=24

2.89e-03

 

N=48

7.91e-04

0.041

5.99e-07

3.74e-08

0.125

1.51e-06

9.45e-08

0.250

2.26e-06

1.41e-07

0.333

2.41e-06

1.51e-07

0.375

2.41e-06

1.50e-07

0.500

2.12e-06

1.32e-07

0.625

1.57e-06

9.86e-08

0.750

9.30e-07

5.81e-08

0.875

3.49e-07

2.18e-08

0.958

8.02e-07

5.01e-09

Table 6. Comparison of the pointwise error $\varepsilon(x)$ when $h=0.05$ and 0.25 for Example 3

$x_i$

Present Method

Reference [13]

Reference [21]

N=10

N=20

N=10

N=20

N=10

8.31e-01

0.5

3.57e-03

2.35e-04

1.28e-01

3.88e-02

1.0

1.77e-03

1.20e-04

3.47e-02

6.56e-02

1.5

8.03e-04

5.70e-05

2.83e-03

9.48e-04

2.0

4.61e-04

3.46e-05

5.30e-03

2.64e-03

2.5

3.57e-04

2.77e-05

7.20e-04

3.03e-04

3.0

3.27e-04

2.57e-05

7.64e-03

3.11e-03

3.5

3.17e-04

2.50e-05

7.73e-03

3.10e-03

4.0

3.03e-04

2.39e-05

7.70e-03

3.01e-03

4.5

2.48e-04

1.95e-05

7.08e-03

2.57e-03

The results are computed as the $l^2(E, h)$, and $l^{\infty}(E, h)$ error norms for the proposed method in Tables 1-3. The error norms of $\varepsilon(x)$ are reported in Tables 4-6 for different space and compared with the results given by references [13, 18, 21, 25, 26] at different space levels $h \leq 1$ are reported for different values of $N$, at domain $[a, b]$. It is clear that from the all tables our results are better than the results given by references [13, 18, 21, 25, 26]. In addition, that the number of subintervals increases, the pointwise errors become smaller. The solution profiles in Figures 1 and 2 visually compare the physical behavior of both the exact and approximate solutions of the problem at different levels of space. Approximate solutions generated through the proposed approach of compact finite difference have a strong agreement with the exact solutions. The consistency across the different levels of space serves to underpin the solver's effectiveness at Fredholm Integro-Differential Equations with high accuracy.

Figure 1. Exact and approximate solution of Example 1 using $N=20, h=0.05$

Figure 2. Exact and approximate solution of Examples 2 and 3 using $N=20, h=0.05$ and $h=0.25$, respectively

Table 7. Rate convergence for Example 4

$N$

$l^2\left(\mathbf{E}_u, h\right)$

Rate

$\boldsymbol{l}^{\infty}\left(\mathbf{E}_u, \boldsymbol{h}\right)$

Rate

$\boldsymbol{l}^2\left(\mathbf{E}_v, \boldsymbol{h}\right)$

Rate

$\boldsymbol{l}^{\infty}\left(\mathbf{E}_{\boldsymbol{v}}, \boldsymbol{h}\right)$

Rate

10

7.2085e-05

 

7.2541e-05

 

1.1759e-04

 

1.0227e-04

 

20

4.5080e-06

3.9991

4.5370e-06

3.9990

7.3236e-06

4.0051

6.3556e-06

4.0082

30

2.8172e-07

4.0002

2.8364e-07

3.9996

4.5719e-07

4.0017

3.9667e-07

4.0020

40

1.7606e-08

4.0001

1.7771e-08

3.9965

2.8565e-08

4.0005

2.4792e-08

4.0000

Figure 3. Exact and approximate solution of u(x) for Example 4 with N=20 and h=0.1571

Figure 4. Exact and approximate solution of v(x) for Example 4 with N=20 and h=0.1571

Table 8. Rate convergence of for Example 5

$N$

$l^2\left(\mathbf{E}_v, h\right)$

Rate

$\boldsymbol{l}^{\infty}\left(\mathbf{E}_u, \boldsymbol{h}\right)$

Rate

$l^2\left(\mathbf{E}_v, h\right)$

Rate

$\boldsymbol{l}^{\infty}\left(\mathbf{E}_v, \boldsymbol{h}\right)$

Rate

10

3.1963e-06

 

4.4469e-06

 

5.0889e-06

 

6.9621e-06

 

20

1.9978e-07

3.9999

2.7793e-07

4.0000

3.1807e-07

3.9999

4.3572e-07

3.9980

40

1.2486e-08

4.0000

1.7411e-08

3.9967

1.9880e-08

4.0000

2.7283e-08

3.9973

80

7.8035e-10

4.0000

1.0881e-09

4.0001

1.2425e-09

4

1.7052e-09

4.0000

Table 9. Comparison between the present method and B-spline for Example 5

$x_i$

Present $\boldsymbol{u}(\boldsymbol{x})$

[19]

Present $v(\boldsymbol{x})$

[19]

N=10

N=10

N=10

N=10

0.1

1.14e-06

3.23e-07

2.73e-06

8.40e-06

0.2

2.22e-06

2.58e-06

4.76e-06

3.70e-05

0.3

3.16e-06

8.73e-06

6.11e-06

9.08e-05

0.4

3.89e-06

2.06e-05

6.83e-06

1.75e-04

0.5

4.34e-06

4.04e-05

6.96e-06

2.95e-04

0.6

4.44e-06

6.98e-05

6.53e-06

4.55e-04

0.7

4.13e-06

1.11e-04

5.57e-06

6.61e-04

0.8

3.33e-06

1.65e-04

4.14e-06

9.17e-04

0.9

1.98e-06

2.35e-04

2.27e-06

1.23e-03

Figure 5. Exact and approximate solution of u(x) for Example 5 with N=20 and h=0.05

Figure 6. Exact and approximate solution of v(x) for Example 5 with N=20 and h=0.05

Table 10. Rate convergence for Example 6

$N$

$l^2\left(\mathrm{E}_u, h\right)$

Rate

$l^{\infty}\left(\mathbf{E}_u, \boldsymbol{h}\right)$

Rate

$\boldsymbol{l}^2\left(\mathbf{E}_v, \boldsymbol{h}\right)$

Rate

$l^{\infty}\left(\mathbf{E}_v, \boldsymbol{h}\right)$

Rate

10

3.5239e-07

 

4.8050e-07

 

2.0983e-07

 

2.8764e-07

 

20

2.2050e-08

3.9983

3.0228e-08

3.9906

1.3125e-08

3.9906

1.8001e-08

3.9981

40

1.3785e-09

3.9996

1.8898e-09

3.9976

8.2048e-10

3.9997

1.1279e-09

3.9964

80

8.6105e-11

4.0009

1.1811e-10

3.9999

5.1280e-11

4

7.0490e-11

4.0001

Figure 7. Exact and approximate solution of u(x) for Example 6 with N=20 and h=0.05

Figure 8. Exact and approximate solution of v(x) for Example 6 with N=20 and h=0.05

To provide the summary of the proposed method to find the approximate solutions based on applying compact finite difference on the system of Examples 4-6 that have been illustrated in Tables 7-10. The error norms of $\mathrm{E}_u(x)$, and $\mathrm{E}_v(x)$ are reported in the Table 9 for space levels and compared with the results given by Ebrahimi and Rashidinia [19]. It is clear that from Table 9 our results are better than the results [19]. One of the reasons is due to the errors produced by the presented scheme are much close to zero and the obtained numerical solutions indicate that the method is reliable and yields results compatible with the previous studies and analytical solutions. In addition, the scheme is shown that is fourth-order convergent in space as evidenced by the results presented in Tables 7, 8, and 10. The solution behavior depicted in Figures 2-8, evidently indicates a comparison between the exact and approximate solutions of the problem at different spatial modes. The approximate solution obtained from the proposed compact finite difference method is highly compatible with the exact solution. The strong agreement for different spatial modes further justifies the efficiency of the method in generating accurate approximate solutions to the coupled Fredholm Integro-Differential Equations.

5. Conclusion

This study establishes a connection between the fourth-order compact finite difference methods with Simpson’s quadrature rule to address Fredholm Integro-Differential Equations (FIDEs) and System of Fredholm Integro-Differential Equations (SFIDEs) of the second order. The precision and efficacy of the proposed scheme are rigorously evaluated through the application to various test problems, with the assessment based on $l^2$ and $l^{\infty}$ error norms at different spatial resolutions. The convergence for the suggegested methods is proved through Sobolev space. Numerical experiments affirm the efficiency, reliability, fruitfulness, and robustness of the presented method in obtaining accurate solutions for SFIDEs. The scheme exhibits a fourth-order spatial accuracy, demonstrating excellent agreement with analytical solutions and outperforming existing solutions reported in the literature.

  References

[1] Cao, H.H., Liu, L.B., Zhang, Y., Fu, S.M. (2011). A fourth-order method of the convection–Diffusion equations with Neumann boundary conditions. Applied Mathematics and Computation, 217(22): 9133-9141. https://doi.org/10.1016/j.amc.2011.03.141

[2] Fu, Y., Tian, Z., Liu, Y. (2018). A compact exponential scheme for solving 1D unsteady convection-diffusion equation with Neumann boundary conditions. arXiv preprint arXiv:1805.05728. https://doi.org/10.48550/arXiv.1805.05728

[3] Yao, C., Zhang, Y., Chen, J., Ling, X., Jing, K., Lu, Y., Fan, E. (2020). Development of a fourth-order compact finite difference scheme for simulation of simulated-moving-bed process. Scientific reports, 10(1): 7820. https://doi.org/10.1038/s41598-020-64562-8

[4] Pirdawood, M.A., Sabawi, Y.A. (2021). High-order solution of Generalized Burgers–Fisher Equation using compact finite difference and DIRK methods. IOP Conference Series: Journal of Physics, 1999(1): 012088. https://doi.org/10.1088/1742-6596/1999/1/012088

[5] Sabawi, Y.A., Pirdawood, M.A., Sadeeq, M.I. (2021). A compact fourth-order implicit-explicit Runge-Kutta type method for solving diffusive Lotka–Volterra system. IOP Conference Series: Journal of Physics, 1999(1): 012103. https://doi.org/10.1088/1742-6596/1999/1/012103

[6] Biazar, J., Mehrlatifan, M.B. (2018). A compact finite difference scheme for reaction-convection-diffusion equation. Chiang Mai Journal of Science, 45(3): 1559-1568. http://www.researchgate.net/publication/316628422.

[7] Biazar, J., Asayesh, R. (2020). An efficient high-order compact finite difference method for the Helmholtz equation. Computational Methods for Differential Equations, 8(3): 553-563. https://doi.org/10.22034/cmde.2020.27993.1382

[8] Roul, P., Goura, V.P., Agarwal, R. (2019). A compact finite difference method for a general class of nonlinear singular boundary value problems with Neumann and Robin boundary conditions. Applied Mathematics and Computation, 350: 283-304. https://doi.org/10.1016/j.amc.2019.01.001

[9] Cai, Y., Fu, J., Liu, J., Wang, T. (2022). A fourth-order compact finite difference scheme for the quantum Zakharov system that perfectly inherits both mass and energy conservation. Applied Numerical Mathematics, 178: 1-24. https://doi.org/10.1016/j.apnum.2022.03.009

[10] Wazwaz, A.M. (2011). Linear and Nonlinear Integral Equations. Springer, Berlin.

[11] Zhou, F., You, Y., Li, G., Xie, G. (2018). The precise integration method for semi-discretized equation in the dual reciprocity method to solve three-dimensional transient heat conduction problems. Engineering Analysis with Boundary Elements, 95: 160-166. https://doi.org/10.1016/j.enganabound.2018.07.005

[12] Sabawi, Y.A., Hussen, O.B. (2024). A cubic B-spline finite element method for Volterra integro-differential equation. Palestine Journal of Mathematics, 13(3): 571-583. 

[13] Jalius, C., Abdul Majid, Z. (2017). Numerical solution of second-order Fredholm integrodifferential equations with boundary conditions by quadrature-difference method. Journal of Applied Mathematics, 2017(1): 2645097. https://doi.org/10.1155/2017/2645097

[14] MacKinnon, R.J., Johnson, R.W. (1991). Differential-equation-based representation of truncation errors for accurate numerical simulation. International Journal for Numerical Methods in Fluids, 13(6): 739-757. https://doi.org/10.1002/fld.1650130606

[15] Soliman, A.F., El-Azab, M.S. (2012). Compact finite difference schemes for partial Integro-Differential Equations. American Academic & Scholarly Research Journal, 4(1): 6-13. https://www.bu.edu.eg/portal/uploads/Engineering,%20Benha/Basic%20Sciences/5821/publications/amal%20foad%20abdelaziz%20soliman_p1.pdf.

[16] Zhao, J., Corless, R.M. (2006). Compact finite difference method for Integro-Differential Equations. Applied mathematics and computation, 177(1): 271-288. https://doi.org/10.1016/j.amc.2005.11.007

[17] Soliman, A.F., El-Asyed, A.M.A., El-Azab, M.S. (2012). Compact finite difference schemes for solving a class of weakly-singular partial Integro-Differential Equations. Mathematical Sciences Letters, 1(1): 53-60. http://www.naturalspublishing.com/files/published/y8x2jzl8k5993m.pdf.

[18] Ebrahimi, N., Rashidinia, J. (2014). Spline collocation for Fredholm and Volterra Integro-Differential Equations. International Journal of Mathematical Modelling & Computations, 4(3 (SUMMER)): 289-298. http://www.sid.ir/FileServer/JE/1033520140305.pdf.

[19] Ebrahimi, N., Rashidinia, J. (2015). Spline collocation for system of Fredholm and Volterra Integro-Differential Equations. Journal of Mathematical Modeling, 3(2): 219-232.

[20] Tair, B., Guebbai, H., Segni, S., Ghiat, M. (2021). Solving linear Fredholm integro-differential equation by Nyström method. Journal of Applied Mathematics and Computational Mechanics, 20(3): 53-64. https://doi.org/10.17512%2Fjamcm.2021.3.05

[21] Hosseini, S.M., Shahmorad, S. (2005). Numerical piecewise approximate solution of Fredholm Integro-Differential Equations by the Tau method. Applied Mathematical Modelling, 29(11): 1005-1021. https://doi.org/10.1016/j.apm.2005.02.003

[22] Elahi, Z., Akram, G., Siddiqi, S.S. (2018). Laguerre approach for solving system of linear Fredholm Integro-Differential Equations. Mathematical Sciences, 12(3): 185-195. https://doi.org/10.1007/s40096-018-0258-0

[23] Seghiri, N., Nadir, M., Khirani, A. (2024). Operational matrices of Genocchi polynomials for solving high-order linear Fredholm Integro-Differential Equations. Mathematical Modelling of Engineering Problems, 11(9): 2473-2482. https://doi.org/10.18280/mmep.110919

[24] Abdulqader, A.J. (2024). Application of the homotopy analysis method in solving fuzzy nonlinear integral equations for birthrate. Mathematical Modelling of Engineering Problems, 11(9): 2527-2536. https://doi.org/10.18280/mmep.110923

[25] Darania, P., Ebadian, A. (2007). A method for the numerical solution of the Integro-Differential Equations. Applied Mathematics and Computation, 188(1): 657-668. https://doi.org/10.1016/j.amc.2006.10.046

[26] Aruchunan, E., Sulaiman, J. (2013). Half-sweep quadrature-difference schemes with iterative method in solving linear Fredholm Integro-Differential Equations. Progress in Applied Mathematics, 5(1): 11-21. https://doi.org/10.3968/j.pam.1925252820130501.2271