Optimal Homotopy Asymptotic and Multistage Optimal Homotopy Asymptotic Methods for Solving System of Volterra Integral Equations of the Second Kind

Abstract

In this paper, optimal homotopy asymptotic method (OHAM) and its implementation on subinterval, called multistage optimal homotopy asymptotic method (MOHAM), are presented for solving linear and nonlinear systems of Volterra integral equations of the second kind. To illustrate these approaches two examples are presented. The results confirm the efficiency and ability of these methods for such equations. The results will be compared to find out which method is more accurate. Advantages of applying MOHAM are also illustrated.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Integral equations, differential equations, integrodifferential equations, and system of such equations, linear and nonlinear, usually appeared in mathematical modeling of different phenomena in physics, biology, and engineering [1–4].

There are various numerical methods for finding approximate solutions of a system of integral equations and a system of integrodifferential equations such as differential transform method [5], Adomian decomposition method [6], modified homotopy perturbation method [7–9], homotopy perturbation method [10], block by block method [11], Galerkin method [12], rationalized Haar functions method [2, 3], Runge-Kutta method [13], homotopy analysis method [14], Tau method [15], and variational iteration method [16].

Moreover, Sezar et al. applied Chebyshev polynomial method and Taylor collocation method for systems of linear differential equations and integrodifferential equations [1, 17]. Yusufoglu employed homotopy pertubation method to solve a system of Fredholm–Volterra type integral equations [18].

One of the powerful and efficient methods for solving integral equations is OHAM. In this paper, we apply OHAM to solve systems of integral equations of the second kind. We will also consider a modified version of OHAM, that is called multistage optimal homotopy asymptotic method (MOHAM). This approach was introduced for the first time by Anakira et al. to approximate the solutions of differential equations with initial-values [19].

The organization of this research is as follows: in Section 2, OHAM and MOHAM are introduced. In Sections 3 and 4, applications of OHAM and MOHAM to system of Volterra integral equations of the second kind are explained, respectively. Section 5 is devoted to proving the convergence of OHAM. In section six, illustrative examples are presented, and conclusion appeared in the last section.

2. OHAM and MOHAM

The OHAM approach is usually applied to solve boundary value problems; say $$\begin{array}{c}\text{(1)}& L_i\left(u_i\left(t\right)\right)+f_i\left(t\right)+N_i\left(u_i\left(t\right)\right)=\mathrm{0},B\left(u_i,\frac{d{u}_i}{dt}\right)=\mathrm{0},\\ \hspace{1em}i=\mathrm{1,2},\dots ,n,\end{array}$$where $L_i$ and $N_i$ are linear and nonlinear functional operators, respectively. $f_i\left(t\right)$ is a known function, $u_i\left(t\right)$ is an unknown function, and $B$ is a boundary operator [20–24].

According to OHAM we construct a homotopy ${\phi }_i\left(u_i\left(t,P\right);P\right)\text{:}\hspace{0.17em}R\times \left[\mathrm{0,1}\right]\to R$ for (1) as follows: $$\begin{array}{c}\text{(2)}& \left(\mathrm{1}-P\right)\left[L_i\left(u_i\left(t,P\right)\right)+f_i\left(t\right)\right]=H\left(P\right)\left[L_i\left(u_i\left(t,P\right)\right)+f_i\left(t\right)+N_i\left(u_i\left(t,P\right)\right)\right],B\left(u_i\left(t,P\right),\frac{d{u}_i\left(t,P\right)}{dt}\right)=\mathrm{0},\end{array}$$where $P\in \left[\mathrm{0,1}\right]$ is an embedding parameter; $H\left(P\right)$, for $P\ne 0$, is a nonzero auxiliary function such that $H\left(0\right)=0$. For $P=0$ and $P=1$ we have $u_i\left(t,0\right)=u_{i0}\left(t\right)$, and $u_i\left(t,1\right)=u_i\left(t\right)$, respectively. Thus, as $P$ increases from $0$ to $1$, the solution $u\left(t,P\right)$ varies from $u_0\left(t\right)$ to the solution $u\left(t\right)$, where $u_0\left(t\right)$ is an initial guess, for the solution, that satisfies the linear operator which is obtained from (2), for $P=0$;$$\begin{array}{c}\text{(3)}& L_i\left(u_{i\mathrm{0}}\left(t\right)\right)+f_i\left(t\right)=\mathrm{0},B\left(u_{i\mathrm{0}},\frac{d{u}_{i\mathrm{0}}}{dt}\right)=\mathrm{0}.\end{array}$$The auxiliary function $H\left(P\right)$ is considered as the following power series in $P$:$$\begin{array}{}\text{(4)}& H\left(P\right)=C_{\mathrm{1}}P+C_{\mathrm{2}}{P}^{\mathrm{2}}+\cdots ,\end{array}$$where $C_1$, $C_2,\dots$ are constants that will be determined later.

Let us consider the approximate solution, $u_i\left(t;P,C_1,C_2,\dots \right)$, as a power series about $P$$$\begin{array}{}\text{(5)}& u_i\left(t;P,C_{\mathrm{1}},C_{\mathrm{2}},\dots \right)=u_{i\mathrm{0}}\left(t\right)+{\displaystyle \sum }_{k\ge \mathrm{1}}{u}_{ik}\left(t,C_{\mathrm{1}},C_{\mathrm{2}},\dots ,C_k\right)P^k,\hspace{1em}i=\mathrm{1,2},\dots ,n.\end{array}$$Substitution in (2) from (5) and equating the coefficients of the terms with identical powers of $P$ lead to governing equations of $u_{i0}\left(t\right)$, $u_{i1}\left(t\right),\dots ,u_{ik}\left(t\right)$, for $i=1,2,\dots ,n$, which starts from (3), followed by$$\begin{array}{c}\text{(6)}& L_i\left(u_{i\mathrm{1}}\left(t\right)\right)=C_{\mathrm{1}}{N}_{i\mathrm{0}}\left(u_{\mathrm{0}}\left(t\right)\right),B\left(u_{i\mathrm{1}},\frac{d{u}_{i\mathrm{1}}}{dt}\right)=\mathrm{0},\\ \text{(7)}& L_i\left(u_{ik}\left(t\right)-u_{ik-\mathrm{1}}\left(t\right)\right)=C_k{N}_{i\mathrm{0}}\left(u_{i\mathrm{0}}\left(t\right)\right)+{\displaystyle \sum }_{l=\mathrm{1}}^{k-\mathrm{1}}{C}_i\left[L_i\left(u_{ik-l}\left(t\right)\right)+N_{ik-l}\left(u_{i\mathrm{0}}\left(t\right),u_{i\mathrm{1}}\left(t\right),\dots ,u_{ik-l}\left(t\right)\right)\right],B\left(u_{ik},\frac{d{u}_{ik}}{dt}\right)=\mathrm{0},\hspace{1em}k=\mathrm{2,3},\dots \end{array}$$where $N_m\left(u_{\mathrm{0}}\left(t\right),u_{\mathrm{1}}\left(t\right),\dots ,u_m\left(t\right)\right)$ is the coefficient of $P^m$ in the expansion of $N\left(u\left(t;P,C_{\mathrm{1}},\dots ,C_m\right)\right)$ about the embedding parameter $P$$$\begin{array}{}\text{(8)}& N_i\left(u_i\left(t;P\right)\right)=N_i\left(u_{i\mathrm{0}}\left(t\right)\right)+{\displaystyle \sum }_{k\ge \mathrm{1}}{N_{ik}\left(u_{i\mathrm{0}}\left(t\right),u_{i\mathrm{1}}\left(t\right),\dots ,u_{ik}\left(t\right)\right)P}^k,\hspace{1em}i=\mathrm{1,2},\dots ,n\end{array}$$where $u_i\left(t;P\right)$ is given by (5).

Studying the rate of convergence of the series (5) depends upon the auxiliary constants $C_q$, $q=\mathrm{1,2},3,\dots$. If the series (5) converges for $P=1$, one has $$\begin{array}{}\text{(9)}& u_i\left(t,C_{\mathrm{1}},C_{\mathrm{2}},\dots \right)=u_{i\mathrm{0}}\left(t\right)+{\displaystyle \sum }_{k\ge \mathrm{1}}{u}_{ik}\left(t,C_{\mathrm{1}},C_{\mathrm{2}},\dots ,C_k\right),\hspace{1em}i=\mathrm{1,2},\dots ,n.\end{array}$$Then the $m$th order approximation is as follows: $$\begin{array}{}\text{(10)}& {\tilde{u}}_i\left(t,C_{\mathrm{1}},C_{\mathrm{2}},\dots ,C_m\right)=u_{i\mathrm{0}}\left(t\right)+{\displaystyle \sum }_{k=\mathrm{1}}^m{u}_{ik}\left(t,C_{\mathrm{1}},C_{\mathrm{2}},\dots ,C_k\right),\hspace{1em}i=\mathrm{1,2},\dots ,n.\end{array}$$Substitution of (10) into (1) results in the following expression for residual: $$\begin{array}{}\text{(11)}& R_i\left(t,C_{\mathrm{1}},C_{\mathrm{2}},\dots ,C_m\right)=L\left({\tilde{u}}_i\left(t,C_{\mathrm{1}},C_{\mathrm{2}},\dots ,C_m\right)\right)+f_i\left(t\right)+N_i\left({\tilde{u}}_i\left(t,C_{\mathrm{1}},C_{\mathrm{2}},\dots ,C_m\right)\right),\hspace{1em}i=\mathrm{1,2},\dots ,n.\end{array}$$If $R_i\left(t,C_1,C_2,\dots ,C_m\right)=0$, then ${\tilde{u}}_i\left(t,C_1,C_2,\dots ,C_m\right)$ will be an exact solution and this, in general, does not happen especially in nonlinear problems. In order to find the optimal values of $C_q$, $q=\mathrm{1,2},3,\dots m$, we apply least square minimization approach $$\begin{array}{}\text{(12)}& \frac{\partial {\mathrm{J}}_{\mathrm{i}}}{\partial {\mathrm{C}}_{\mathrm{1}}}=\frac{\partial {\mathrm{J}}_{\mathrm{i}}}{\partial {\mathrm{C}}_{\mathrm{2}}}=\cdots =\frac{\partial {\mathrm{J}}_{\mathrm{i}}}{\partial {\mathrm{C}}_{\mathrm{m}}}=\mathrm{0},\end{array}$$where$$\begin{array}{}\text{(13)}& J_i\left(C_{\mathrm{1}},C_{\mathrm{2}},\dots ,C_m\right)={{\displaystyle \int }}_a^b{R}_i^{\mathrm{2}}\left(t,C_{\mathrm{1}},C_{\mathrm{2}},\dots ,C_m\right)dt,\hspace{1em}i=\mathrm{1,2},\dots ,n,\end{array}$$where $a$ and $b$ are two values, depending on the given problem. Knowing $C_q$, $q=\mathrm{1,2},\mathrm{3},\dots m$, from (12), the approximate solution of order $m$ will be determined easily.

If the interval of changes of the time variable is long, then OHAM fails to reach accurate solutions.

MOHAM overcomes this shortcoming by partitioning the time interval, $\left[t_0,T\right]$, into $N$ subintervals $\left[t_0,t_1\right),\dots ,\left[t_{\gamma -1},t_{\gamma }\right]$, where $t_{\gamma }=T$ and OHAM will be applied over each subintervals. The solution at the last point, in each subinterval, denotes an initial approximation to the solution, over the next interval. The process will continue until we achieve the preassigned time, $T$

Implementation of MOHAM is almost the same as OHAM, with some minor changes.

Equations (4), (10), (11), (12), and (13) change to (16), (17), (18), (20), and (19), respectively. Also, initial approximation in ${[t}_{\gamma },t_{\gamma +1})$, $\gamma =0,1,\dots ,N-1$, will be considered as $$\begin{array}{}\text{(14)}& u_{i\mathrm{0}}^j\left(t_j\right)={\alpha }_i^j,\hspace{1em}i=\mathrm{1,2},\dots ,n,\hspace{0.17em}\hspace{0.17em}j=\mathrm{1},\dots ,N,\end{array}$$In addition, deformation equation in each subinterval for $i=\mathrm{1,2},\dots ,n$, $j=1,\dots ,N$, will change to the following [19]:$$\begin{array}{}\text{(15)}& \left(\mathrm{1}-P\right)\left[L_i\left(u_i^j\left(t,P\right)\right)-u_{i\mathrm{0}}^j\left(t\right)\right]=H\left(P\right)\left[L_i\left(u_i^j\left(t,P\right)\right)+f_i\left(t\right)+N_i\left(u_i^j\left(t,P\right)\right)\right].\end{array}$$Moreover, auxiliary function $H^j\left(P,t\right)$ will be generalized as follows:$$\begin{array}{}\text{(16)}& H^j\left(P,t\right)=\left(C_{i\mathrm{1}}^j+C_{i\mathrm{2}}^{j}t+C_{i\mathrm{3}}^j{t}^{\mathrm{2}}+\cdots \right)P.\text{(17)}& {\tilde{u}}_i^j\left(t,C_{iq}^j\right)=u_{i\mathrm{0}}^j\left(t\right)+{\displaystyle \sum }_{k=\mathrm{1}}^m{u}_{k,j}\left(t,C_{iq}^j\right),\hspace{1em}q=\mathrm{1,2},\mathrm{3},\dots ,m,\text{(18)}& R_i^j\left(t,C_{iq}^j\right)=L\left({\tilde{u}}_i\left(t,C_{iq}^j\right)\right)+f_i\left(t\right)+N_i\left({\tilde{u}}_i\left(t,C_{iq}^j\right)\right),\hspace{1em}q=\mathrm{1,2},\mathrm{3},\dots ,m,\text{(19)}& J_i^j\left(C_{iq}^j\right)={{\displaystyle \int }}_{t_{\gamma }}^{t_{\gamma }+h}{\left(R_i^j\right)}^{\mathrm{2}}\left(s,C_{iq}^j\right)ds,\hspace{1em}\gamma =\mathrm{0},\mathrm{1},\dots ,N-\mathrm{1},\hspace{0.17em}\hspace{0.17em}q=\mathrm{1,2},\mathrm{3},\dots ,m.\end{array}$$The length of the subinterval $\left[t_{\gamma },t_{\gamma +1}\right)$ is apparently $h$, and the number of subintervals is $N=⌊T/h⌋$. Now, we consider derivatives of (19), with respect to $C_{iq}^j$, $(i=\mathrm{1,2},\dots ,n,j=1,\dots ,N,q=\mathrm{1,2},3,\dots ,m)$ to zero. In fact we define ${\alpha }_i^j={\tilde{u}}_i^j\left(t_j\right)$, in each subinterval $\left[t_{\gamma },t_{\gamma +1}\right)$. Therefore, the convergence control parameters can be determined from the solution of the following system of equations:$$\begin{array}{}\text{(20)}& \frac{\partial J^j}{\partial {\mathrm{C}}_{\mathrm{i}\mathrm{1}}^{\mathrm{j}}}=\frac{\partial J^j}{\partial {\mathrm{C}}_{\mathrm{i}\mathrm{2}}^{\mathrm{j}}}=\cdots =\frac{\partial J^j}{\partial {\mathrm{C}}_{\mathrm{i}\mathrm{m}}^{\mathrm{j}}}=\mathrm{0}.\end{array}$$Approximate analytic solutions, on each subinterval, are as follows:$$\begin{array}{}\text{(21)}& \tilde{u}\left(t\right)=\left\{\begin{array}{cc}{\tilde{u}}_{\mathrm{1}}\left(t\right),\hfill & t_{\mathrm{0}}\le t<t_{\mathrm{1}},\hfill \\ {\tilde{u}}_{\mathrm{2}}\left(t\right),\hfill & t_{\mathrm{1}}\le t<t_{\mathrm{2}},\hfill \\ ⋮\hfill & \hfill \\ {\tilde{u}}_N\left(t\right),\hfill & t_{N-\mathrm{1}}\le t\le T.\hfill \end{array}\right.\end{array}$$

3. Application of OHAM to Systems of Volterra Integral Equations

In this section, we apply OHAM on the following system of Volterra integral equations: $$\begin{array}{}\text{(22)}& U\left(t\right)=F\left(t\right)+{{\displaystyle \int }}_{\mathrm{0}}^{t}K\left(t,s\right)U\left(s\right)ds,\hspace{1em}t\in \left[\mathrm{0,1}\right],\end{array}$$where $t\in R$ and$$\begin{array}{c}\text{(23)}& K\left(t,s\right)=\left[k_{i,j}\left(t,s\right)\right],\hspace{1em}i,j=\mathrm{1,2},\dots ,n,U\left(t\right)={\left[u_{\mathrm{1}}\left(t\right),u_{\mathrm{2}}\left(t\right),\dots ,u_n\left(t\right)\right]}^t,\\ F\left(t\right)={\left[f_{\mathrm{1}}\left(t\right),f_{\mathrm{2}}\left(t\right),\dots ,f_n\left(t\right)\right]}^t.\end{array}$$In (22), $K$ is the kernel that is known, $F$ is given vector, and $U$ is an unknown vector function. Consider $i$ th equation of (22)$$\begin{array}{}\text{(24)}& u_i\left(t\right)=f_i\left(t\right)+{{\displaystyle \int }}_{\mathrm{0}}^t{\displaystyle \sum }_{j=\mathrm{1}}^n{k}_{i,j}\left(t,s\right)u_j\left(s\right)ds,\hspace{1em}i=\mathrm{1,2},\cdots ,n.\end{array}$$Using aforementioned OHAM procedure results in the following sequential equations, for $i=\mathrm{1,2},\dots ,n.$$$\begin{array}{c}\text{(25)}& P^{\mathrm{0}}\text{:}{\hspace{0.17em}\hspace{0.17em}u}_{i\mathrm{0}}\left(t\right)=f_i\left(t\right),P^{\mathrm{1}}\text{:}{\hspace{0.17em}\hspace{0.17em}u}_{i\mathrm{1}}\left(t\right)=-C_{i\mathrm{1}}{{\displaystyle \int }}_{\mathrm{0}}^t{\displaystyle \sum }_{j=\mathrm{1}}^n{k}_{i,j}\left(t,s\right)u_{j\mathrm{0}}\left(s\right)ds,\\ P^{\mathrm{2}}\text{:}{\hspace{0.17em}\hspace{0.17em}u}_{i\mathrm{2}}\left(t\right)=\left(\mathrm{1}+C_{i\mathrm{1}}\right)u_{i\mathrm{1}}\left(t\right)-C_{i\mathrm{1}}{{\displaystyle \int }}_{\mathrm{0}}^t{\displaystyle \sum }_{j=\mathrm{1}}^n{k}_{i,j}\left(t,s\right)u_{j\mathrm{1}}\left(s\right)ds-C_{i\mathrm{2}}{{\displaystyle \int }}_{\mathrm{0}}^t{\displaystyle \sum }_{j=\mathrm{1}}^n{k}_{i,j}\left(t,s\right)u_{j\mathrm{0}}\left(s\right)ds,\\ ⋮\\ P^k\text{:}{\hspace{0.17em}\hspace{0.17em}u}_{ik}\left(t\right)=\left(\mathrm{1}+C_{i\mathrm{1}}\right)u_{ik-l}\left(t\right)+{\displaystyle \sum }_{l=\mathrm{2}}^{k-\mathrm{1}}{C}_{il}{u}_{ik-l}\left(t\right)-{\displaystyle \sum }_{z=\mathrm{1}}^k{C}_{iz}{{\displaystyle \int }}_{\mathrm{0}}^t{\displaystyle \sum }_{j=\mathrm{1}}^n{k}_{i,j}\left(t,s\right)u_{jk-z}\left(s\right)ds,\\ \hspace{1em}k=\mathrm{3,4},\dots .\end{array}$$Using (11)–(13), we find $C_{\beta }$, $\beta =\mathrm{1},\dots ,m.$

Knowing these parameters, an approximate solution, of order $m$, will be determined.

4. Application of MOHAM to Systems of Volterra Integral Equations

In this section, we apply MOHAM to (22). This procedure leads to the following sequence of equations, for $i=\mathrm{1,2},\dots ,n.$$$\begin{array}{c}\text{(26)}& P^{\mathrm{0}}\text{:}{\hspace{0.17em}\hspace{0.17em}u}_{i\mathrm{0}}^j\left(t\right)={\alpha }_i^j,P^{\mathrm{1}}\text{:}{\hspace{0.17em}\hspace{0.17em}u}_{i\mathrm{1}}^j\left(t\right)=\left(C_{i\mathrm{1}}^j+C_{i\mathrm{2}}^{j}t+\cdots \right)\left(u_{i\mathrm{0}}^j\left(t\right)-f_i\left(t\right)-{{\displaystyle \int }}_{\mathrm{0}}^t{\displaystyle \sum }_{j=\mathrm{1}}^n{k}_{i,j}{u}_{j\mathrm{0}}^j\left(t\right)ds\right),\\ P^{\mathrm{2}}\text{:}{\hspace{0.17em}\hspace{0.17em}u}_{i\mathrm{2}}^j\left(t\right)=u_{i\mathrm{1}}^j\left(t\right)+\left(C_{i\mathrm{1}}^j+C_{i\mathrm{2}}^{j}t+\cdots \right)\left(u_{i\mathrm{1}}^j\left(t\right)-{{\displaystyle \int }}_{\mathrm{0}}^t{\displaystyle \sum }_{j=\mathrm{1}}^n{k}_{i,j}\left(t,s\right)u_{j\mathrm{1}}^j\left(t\right)ds\right),\\ ⋮\\ P^k\text{:}{\hspace{0.17em}\hspace{0.17em}u}_{j\mathrm{1}}^j\left(t\right)=u_{ik-\mathrm{1}}^j\left(t\right)+\left({C^j}_{i\mathrm{1}}+{C^j}_{i\mathrm{2}}t+\cdots \right)\left(u_{ik-\mathrm{1}}^j\left(t\right)-{{\displaystyle \int }}_{\mathrm{0}}^t{\displaystyle \sum }_{j=\mathrm{1}}^n{k}_{i,j}\left(t,s\right)u_{jk-\mathrm{1}}^j\left(t\right)ds\right),\\ \hspace{1em}k=\mathrm{3,4},\dots .\end{array}$$In addition, by using (18)–(20), we find $C_{iq}^j$, $i=\mathrm{1,2},\mathrm{3},\dots n$, $j=\mathrm{1},\dots ,N$, $q=\mathrm{1,2},\mathrm{3},\dots ,m.$

Knowing that the values are parameters, the approximate solution of order $m$ will be determined.

5. Convergence of the OHAM

There are two proofs presented in [19, 20] that have some oversight. The following proof covers these shortcoming.

Theorem 1.

If series (9) convergences to $u_i(t)$ where $u_{ik}\left(t,C_1,\dots ,C_k\right)\in L^2\left[c,d\right]$ is produced by (3), (6), and (7), then $u_i(t)$ is the exact solution of (1).

Proof.

Since the series$$\begin{array}{}\text{(27)}& {\displaystyle \sum }_{k=\mathrm{0}}^{\infty }{u}_{ik}\left(t,C_{\mathrm{1}},\dots ,C_k\right)\end{array}$$converges; it holds$$\begin{array}{}\text{(28)}& \underset{k\to \infty }{\lim }{u}_{ik}\left(t,C_{\mathrm{1}},\dots ,C_k\right)=\mathrm{0}.\end{array}$$The left hand-side of (7) satisfies$$\begin{array}{}\text{(29)}& u_{i\mathrm{1}}\left(t,C_{\mathrm{1}}\right)+{\displaystyle \sum }_{k=\mathrm{2}}^n{u}_{ik}\left(t,C_{\mathrm{1}},\dots ,C_k\right)-{\displaystyle \sum }_{k=\mathrm{2}}^n{u}_{ik-\mathrm{1}}\left(t,{C_{\mathrm{1}},\dots ,C}_{k-\mathrm{1}}\right){=u}_{in}\left(t,C_{\mathrm{1}},\dots ,C_n\right).\end{array}$$According to (29), we have $$\begin{array}{}\text{(30)}& u_{i\mathrm{1}}\left(t,C_{\mathrm{1}}\right)+{\displaystyle \sum }_{k=\mathrm{2}}^{\infty }{u}_{ik}\left(t,C_{\mathrm{1}},\dots ,C_k\right)-{\displaystyle \sum }_{k=\mathrm{2}}^{\infty }{u}_{ik-\mathrm{1}}\left(t,{C_{\mathrm{1}},\dots ,C}_{k-\mathrm{1}}\right)=\underset{n\to \infty }{\lim }{u}_{in}\left(t,C_{\mathrm{1}},\dots ,C_n\right)=\mathrm{0}.\end{array}$$Since the operator $L_i$ is linear, from (30), we have $$\begin{array}{}\text{(31)}& L_i\left(u_{i\mathrm{1}}\left(t,C_{\mathrm{1}}\right)\right)+{\displaystyle \sum }_{k=\mathrm{2}}^{\infty }{L}_i\left(u_{ik}\left(t,C_{\mathrm{1}},\dots ,C_k\right)\right)-{\displaystyle \sum }_{k=\mathrm{2}}^{\infty }{L}_i\left(u_{ik-\mathrm{1}}\left(t,{C_{\mathrm{1}},\dots ,C}_{k-\mathrm{1}}\right)\right)=L_i\left(u_{i\mathrm{1}}\left(t,C_{\mathrm{1}}\right)\right)+L_i\left({\displaystyle \sum }_{k=\mathrm{2}}^{\infty }{u}_{ik}\left(t,C_{\mathrm{1}},\dots ,C_k\right)\right)-L_i\left({\displaystyle \sum }_{k=\mathrm{2}}^{\infty }{u}_{ik-\mathrm{1}}\left(t,{C_{\mathrm{1}},\dots ,C}_{k-\mathrm{1}}\right)\right)=\mathrm{0}.\end{array}$$Then$$\begin{array}{}\text{(32)}& L_i\left(u_{i\mathrm{1}}\left(t,C_{\mathrm{1}}\right)\right)+L_i\left({\displaystyle \sum }_{k=\mathrm{2}}^{\infty }{u}_{ik}\left(t,C_{\mathrm{1}},\dots ,C_k\right)\right)-L_i\left({\displaystyle \sum }_{k=\mathrm{2}}^{\infty }{u}_{ik-\mathrm{1}}\left(t,{C_{\mathrm{1}},\dots ,C}_{k-\mathrm{1}}\right)\right)=L_i\left(u_{i\mathrm{1}}\left(t,C_{\mathrm{1}}\right)\right)+{\displaystyle \sum }_{k=\mathrm{2}}^{\infty }\left(C_k{N}_{i\mathrm{0}}\left(u_{i\mathrm{0}}\left(t\right)\right)+{\displaystyle \sum }_{m=\mathrm{1}}^{k-\mathrm{1}}{C}_m\left[L_i\left(u_{ik-m}\left(t\right)\right)+N_{ik-m}\left(u_{i\mathrm{0}}\left(t\right),u_{i\mathrm{1}}\left(t\right),\dots ,u_{ik-m}\left(t\right)\right)\right]\right)=\mathrm{0},\end{array}$$Equation (32) can be written as follows:$$\begin{array}{}\text{(33)}& L_i\left(u_{i\mathrm{0}}\left(t\right)\right)+f_i\left(t\right)+C_{\mathrm{1}}{N}_{i\mathrm{0}}\left(u_{i\mathrm{0}}\left(t\right)\right)+{\displaystyle \sum }_{k=\mathrm{2}}^{\infty }\left(C_k{N}_{i\mathrm{0}}\left(u_{i\mathrm{0}}\left(t\right)\right)\right.+\left.{\displaystyle \sum }_{m=\mathrm{1}}^{k-\mathrm{1}}{C}_m\left[L_i\left(u_{ik-m}\left(t\right)\right)+N_{ik-m}\left(u_{i\mathrm{0}}\left(t\right),u_{i\mathrm{1}}\left(t\right),\dots ,u_{ik-m}\left(t\right)\right)\right]\right)={\displaystyle \sum }_{k=1}^{\infty }\left({\displaystyle \sum }_{m=\mathrm{1}}^k{C}_{k-m}\left[L\left(u_{m-\mathrm{1}}\left(t,C_{\mathrm{1}},\dots ,C_{m-\mathrm{1}}\right)\right)+N_{im-\mathrm{1}}\left(u_{i\mathrm{0}}\left(t\right),u_{i\mathrm{1}}\left(t\right),\dots ,u_{ik-m}\left(t\right)\right)\right]\right)+f_i\left(t\right)=\mathrm{0}.\end{array}$$If the $C_n$, $n=\mathrm{0,1},\dots ,m$, chosen properly, then (33) leads to (1).

6. Numerical Examples

In this section, two systems of Volterra integral equations of the second kind, a linear and a nonlinear, will be solved to show the efficiency of both OHAM and MOHAM. The results of applying OHAM and MOHAM will be compared. Matlab package is used to carry out computations, with double precision.

Example 1.

let us consider the following linear system of Volterra integral equations of the second kind, with exact solutions $u_1\left(t\right)=e^t$ and $u_2\left(t\right)=e^{-t}$.$$\begin{array}{c}\text{(34)}& u_{\mathrm{1}}\left(t\right)=\mathrm{1}-\frac{t^{\mathrm{2}}}{\mathrm{2}}+{{\displaystyle \int }}_{\mathrm{0}}^t\left(u_{\mathrm{1}}\left(s\right)+s{e}^s{u}_{\mathrm{2}}\left(s\right)\right)ds,u_{\mathrm{2}}\left(t\right)=\mathrm{1}+\frac{t^{\mathrm{2}}}{\mathrm{2}}-{{\displaystyle \int }}_{\mathrm{0}}^t\left(s{e}^{-s}{u}_{\mathrm{1}}\left(s\right)+u_{\mathrm{2}}\left(s\right)\right)ds.\end{array}$$The following aforementioned OHAM procedure results in(35) u 10 t = - t 2 2 + 1 , u 20 t = t 2 2 + 1 , u 11 t = - C 11 t + 0.5 e t t 3 - 3 t 2 + 6 t - 6 + e t t - 1 - 0.1666666667 t 3 + 4 ) , u 21 t = - C 21 2 e - t t + 1 - 2 , u 12 t = 0.041666666 6 C 11 2 - t 4 + 12 t 3 e t + 16 t 3 - 72 t 2 e t + 36 t 2 + 192 t e t + 96 t - 288 e t + 288 - C 12 + C 11 C 11 + 1 t + 0.5 e t t 3 - 3 t 2 + 6 t - 6 + e t t - 1 - 0.1666666667 t 3 + 4 ) , u 22 t = - C 21 2 t 5 + 0.1666666667 e - t - 0.375 t 4 + 0.6666666667 e - t + 1.333333333 t 3 - e - 2 t + e - t - 2 t 2 + - 2 t - 1 e - 2 t + 1 - C 22 2 e - t t + 1 - 2 - C 22 C 22 + 1 2 e t t + 1 - 2 , u 13 t = C 11 C 11 2 0.02777777778 + 0.1 e t t 6 + C 11 2 0.1416666667 - 0.975 e t t 5 + 6.208333333 C 1 2 e t + 0.04166666667 C 11 C 11 - 1 - C 12 t 4 + 0.1666666667 C 11 C 11 + 1 + C 12 + 0.5 e t C 11 + C 12 + C 11 2 - 26.83333333 e t + e - t t 3 + C 11 2 - 0.5 + 82 e t + 5 e - t + 3 C 11 + C 12 0.5 - e t t 2 + 8 e t C 11 + C 12 + C 11 2 11 e - t - 169 e t + 4 C 11 + C 12 - 8 C 11 2 t + C 11 2 177 e t + 11 e - t - 188 + 12 C 11 + C 12 1 - e t + C 12 - C 12 C 11 + 1 - C 13 - C 11 C 12 t + 0.5 e t t 3 - 3 t 2 + 6 t - 6 + e t t - 1 - 0.1666666667 t 3 + 4 + 0.04166666667 C 11 C 11 C 11 + 1 + C 12 - t 4 + 4 3 e t + 4 t 3 + 36 t 2 1 - 2 e t + 192 t + 96 t 1 + 2 e t - 288 e t + 288 , (36) u 23 t = C 23 t - e - t t + 1 + 0.5 e - t t 3 + 3 t 2 + 6 t + 6 + 0.1666666667 t 3 - 2 - C 21 + 1 C 21 C 21 t 5 - 0.375 C 21 t 4 + 1.333333333 C 21 t 3 - 2 C 21 t 2 + C 21 - 2 C 21 t e - 2 t - C 21 e - 2 t + 0.1666666667 C 21 t 4 e - t + 0.6666666667 C 21 t 3 e - t + C 21 t 2 e - t - e - 2 t + C 22 2 e - t t + 1 - 2 + C 21 C 21 + 1 2 e - t . t + 1 - 2 - C 21 ( 0.01666666667 C 21 2 t 6 - C 21 2 0.04166666667 e - t + 0.075 t 5 + C 21 + C 22 t 5 + C 21 2 0.4583333333 e - t + 0.7083333333 t 4 + C 21 2 24 e - t + 1.25 e - 2 t - 25.25 - 2 C 21 + C 22 - 2 C 21 2 t 2 - C 21 2 t + C 21 + C 22 0.1666666667 e - t - 0.375 t 4 + 0.6666666667 e - t + 1.333333333 t 3 + t 2 e - t - 4 t e - t - 6 e - t - 2 t + 6 + 26 C 21 2 t e - t + 1.5 C 21 2 t e - 2 t + C 21 2 2.666666667 e - t - 2 t 3 + 10 e - t + 0.5 e - 2 t t 2 - C 21 C 22 2 e - t t + 1 - 2 - 0.008333333333 C 21 C 22 e - t e t 12 t 5 - 45 t 4 + 160 t 3 - 240 t 2 - 240 t + 720 + 20 t 4 + 80 t 3 + 120 t 2 - 480 t - 720 ,Replacement of the first, second, third, and forth terms in $u_i^3={\sum }_{m=0}^3{u}_{im}$, $i=\mathrm{1,2}$, results in $$\begin{array}{}\text{(37)}& u_{\mathrm{1}}^{\mathrm{3}}\left(t\right)=C_{\mathrm{11}}\left(C_{\mathrm{11}}^{\mathrm{2}}\left(\mathrm{0.02777777778}+\mathrm{0.1}{e}^t\right){t^{\mathrm{6}}+C}_{\mathrm{11}}^{\mathrm{2}}\left(\mathrm{0.1416666667}-\mathrm{0.975}{e}^t\right)t^{\mathrm{5}}+\left(\mathrm{6.208333333}{C}_{\mathrm{1}}^{\mathrm{2}}{e}^t+\mathrm{0.04166666667}\left(C_{\mathrm{11}}\left(C_{\mathrm{11}}-\mathrm{1}\right)-C_{\mathrm{12}}\right)\right)t^{\mathrm{4}}+\left(\mathrm{0.1666666667}\left(C_{\mathrm{11}}\left(C_{\mathrm{11}}+\mathrm{1}\right)+C_{\mathrm{12}}\right)+\mathrm{0.5}{e}^t\left(C_{\mathrm{11}}+C_{\mathrm{12}}\right)+C_{\mathrm{11}}^{\mathrm{2}}\left(-\mathrm{26.83333333}{e}^t+e^{-t}\right)\right)t^{\mathrm{3}}+\left(C_{\mathrm{11}}^{\mathrm{2}}\left(-\mathrm{0.5}+\mathrm{82}{e}^t+\mathrm{5}{e}^{-t}\right)+\mathrm{3}\left(C_{\mathrm{11}}+C_{\mathrm{12}}\right)\left(\mathrm{0.5}-e^t\right)\right)t^{\mathrm{2}}+\left(\mathrm{8}{e}^t\left(C_{\mathrm{11}}+C_{\mathrm{12}}\right)+C_{\mathrm{11}}^{\mathrm{2}}\left(\mathrm{11}{e}^{-t}-\mathrm{169}{e}^t\right)+\mathrm{4}\left(C_{\mathrm{11}}+C_{\mathrm{12}}\right)-\mathrm{8}{C}_{\mathrm{11}}^{\mathrm{2}}\right)t+C_{\mathrm{11}}^{\mathrm{2}}\left(\left(\mathrm{177}{e}^t+\mathrm{11}{e}^{-t}\right)-\mathrm{188}\right)+\mathrm{12}\left(C_{\mathrm{11}}+C_{\mathrm{12}}\right)\left(\mathrm{1}-e^t\right)\right)+\left(\left(C_{\mathrm{12}}-C_{\mathrm{12}}\right)\left(C_{\mathrm{11}}+\mathrm{1}\right)-C_{\mathrm{13}}-C_{\mathrm{11}}{C}_{\mathrm{12}}\right)\left(t+\mathrm{0.5}{e}^t\left(t^{\mathrm{3}}-\mathrm{3}{t}^{\mathrm{2}}+\mathrm{6}t-\mathrm{6}\right)+e^t\left(t-\mathrm{1}\right)-\mathrm{0.1666666667}{t}^{\mathrm{3}}+\mathrm{4}\right)+\mathrm{0.04166666667}{C}_{\mathrm{11}}\left(C_{\mathrm{11}}\left(C_{\mathrm{11}}+\mathrm{1}\right)+C_{\mathrm{12}}\right)\left(-t^{\mathrm{4}}+\mathrm{4}\left({\mathrm{3}e}^t+\mathrm{4}\right)t^{\mathrm{3}}+\mathrm{36}{t}^{\mathrm{2}}\left(\mathrm{1}-\mathrm{2}{e}^t\right)+\mathrm{192}t+\mathrm{96}t\left(\mathrm{1}+\mathrm{2}{e}^t\right)-\mathrm{288}{e}^t+\mathrm{288}\right)+\mathrm{0.041666666}{\mathrm{6}C}_{\mathrm{11}}^{\mathrm{2}}\left(-t^{\mathrm{4}}+\mathrm{12}{t}^{\mathrm{3}}{e}^t+\mathrm{16}{t}^{\mathrm{3}}-\mathrm{72}{t}^{\mathrm{2}}{e}^t+\mathrm{36}{t}^{\mathrm{2}}+\mathrm{192}t{e}^t+\mathrm{96}t-\mathrm{288}{e}^t+\mathrm{288}\right)-\left(C_{\mathrm{12}}+C_{\mathrm{11}}\left(C_{\mathrm{11}}+\mathrm{1}\right)\right)\left(t+\mathrm{0.5}{e}^t\right)\left(t^{\mathrm{3}}-\mathrm{3}{t}^{\mathrm{2}}+\mathrm{6}t-\mathrm{6}\right)+e^t\left(t-\mathrm{1}\right)-\mathrm{0.1666666667}{t}^{\mathrm{3}}+\mathrm{4}\text{)}-C_{\mathrm{11}}\left(t+\mathrm{0.5}{e}^t\right)\left(t^{\mathrm{3}}-\mathrm{3}{t}^{\mathrm{2}}+\mathrm{6}t-\mathrm{6}\right)+e^t\left(t-\mathrm{1}\right)-\mathrm{0.1666666667}{t}^{\mathrm{3}}+\mathrm{4}\text{)}-\mathrm{0.5}{t}^{\mathrm{2}}+\mathrm{1},\text{(38)}& u_{\mathrm{2}}^{\mathrm{3}}\left(t\right)=C_{\mathrm{23}}\left(t-e^{-t}\left(t+\mathrm{1}\right)+\mathrm{0.5}{e}^{-t}\left(t^{\mathrm{3}}+\mathrm{3}{t}^{\mathrm{2}}+\mathrm{6}t+\mathrm{6}\right)+\mathrm{0.1666666667}{t}^{\mathrm{3}}-\mathrm{2}\right)-\left(C_{\mathrm{21}}+\mathrm{1}\right)\left(C_{\mathrm{21}}\left({C_{\mathrm{21}}t}^{\mathrm{5}}-\mathrm{0.375}{C_{\mathrm{21}}t}^{\mathrm{4}}+\mathrm{1.333333333}{C_{\mathrm{21}}t}^{\mathrm{3}}-\mathrm{2}{C}_{\mathrm{21}}{t}^{\mathrm{2}}+C_{21}-\mathrm{2}{C}_{\mathrm{21}}t{e}^{-\mathrm{2}t}-C_{\mathrm{21}}{e}^{-\mathrm{2}t}\right)+\mathrm{0.1666666667}{C}_{\mathrm{21}}{t}^{\mathrm{4}}{e}^{-t}+\mathrm{0.6666666667}{C}_{\mathrm{21}}{t}^{\mathrm{3}}{e}^{-t}+C_{\mathrm{21}}{t}^{\mathrm{2}}\left(e^{-t}-e^{-\mathrm{2}t}\right)\right)+C_{\mathrm{22}}\left(\mathrm{2}{e}^{-t}\left(t+\mathrm{1}\right)-\mathrm{2}\right)+C_{\mathrm{21}}\left(C_{\mathrm{21}}+\mathrm{1}\right)\left(\mathrm{2}{e}^{-t}\right).\left(t+\mathrm{1}\right)-\mathrm{2}\text{))}-C_{\mathrm{21}}\text{(}\mathrm{0.01666666667}{C}_{\mathrm{21}}^{\mathrm{2}}{t}^{\mathrm{6}}-C_{\mathrm{21}}^{\mathrm{2}}\left(\mathrm{0.04166666667}{e}^{-t}+\mathrm{0.075}\right)t^{\mathrm{5}}+\left(C_{\mathrm{21}}+C_{\mathrm{22}}\right)t^{\mathrm{5}}+C_{\mathrm{21}}^{\mathrm{2}}\left(\mathrm{0.4583333333}{e}^{-t}+\mathrm{0.7083333333}\right)t^{\mathrm{4}}+C_{\mathrm{21}}^{\mathrm{2}}\left(\mathrm{24}{e}^{-t}+\mathrm{1.25}{e}^{-\mathrm{2}t}-\mathrm{25.25}\right)-\mathrm{2}\left(C_{\mathrm{21}}{+C}_{\mathrm{22}}-\mathrm{2}{C}_{\mathrm{21}}^{\mathrm{2}}\right)t^{\mathrm{2}}-C_{\mathrm{21}}^{\mathrm{2}}t+\left(C_{\mathrm{21}}+C_{\mathrm{22}}\right)\left(\left(\mathrm{0.1666666667}{e}^{-t}-\mathrm{0.375}\right)t^{\mathrm{4}}+\left(\mathrm{0.6666666667}{e}^{-t}+\mathrm{1.333333333}\right)t^{\mathrm{3}}+t^{\mathrm{2}}{e}^{-t}-\mathrm{4}t{e}^{-t}-\mathrm{6}{e}^{-t}-\mathrm{2}t+\mathrm{6}\right)+\mathrm{26}{C}_{\mathrm{21}}^{\mathrm{2}}t{e}^{-t}+\mathrm{1.5}{C}_{\mathrm{21}}^{\mathrm{2}}t{e}^{-\mathrm{2}t}+C_{\mathrm{21}}^{\mathrm{2}}\left(\left(\mathrm{2}{e}^{-t}-\mathrm{3.333333333}\right)t^{\mathrm{3}}+\left(\mathrm{9}{e}^{-t}-\mathrm{0.5}{e}^{-\mathrm{2}t}+\mathrm{2}\right)t^{\mathrm{2}}\right)-C_{\mathrm{21}}{C}_{\mathrm{22}}\left(\mathrm{2}{e}^{-t}\left(t+\mathrm{1}\right)-\mathrm{2}\right)-\mathrm{0.008333333333}{C}_{\mathrm{21}}{C}_{\mathrm{22}}{e}^{-t}\left(e^t\left(\mathrm{12}{t}^{\mathrm{5}}-\mathrm{45}{t}^{\mathrm{4}}+\mathrm{160}{t}^{\mathrm{3}}-\mathrm{240}{t}^{\mathrm{2}}-\mathrm{240}t+\mathrm{720}\right)+\mathrm{20}{t}^{\mathrm{4}}+\mathrm{80}{t}^{\mathrm{3}}+\mathrm{120}{t}^{\mathrm{2}}-\mathrm{480}t-\mathrm{720}\right)-C_{\mathrm{21}}^{\mathrm{2}}\left(t^{\mathrm{5}}+\left(\mathrm{0.1666666667}{e}^{-t}-\mathrm{0.375}\right)t^{\mathrm{4}}+\left(-\mathrm{2}t-\mathrm{1}\right)e^{-\mathrm{2}t}+\mathrm{1}\right)-C_{\mathrm{22}}\left(\left(C_{\mathrm{22}}+\mathrm{1}\right)\left(\mathrm{2}{e}^t\left(t+\mathrm{1}\right)-\mathrm{2}\right)+\mathrm{2}{e}^{-t}\left(t+\mathrm{1}\right)-\mathrm{2}\right)-C_{\mathrm{21}}\left(\mathrm{2}{e}^{-t}\left(t+\mathrm{1}\right)-\mathrm{2}\right)+\mathrm{0.5}{t}^{\mathrm{2}}+\mathrm{1}.\end{array}$$Using the same technique as presented by (11)–(13), we find $C_{i1}$, $C_{i2}$, and $C_{i3}$, $i=\mathrm{1,2}$, as follows:$$\begin{array}{c}\text{(39)}& C_{\mathrm{11}}=-\mathrm{1.1515281477},C_{\mathrm{12}}=\mathrm{0.0312808541},\\ C_{\mathrm{13}}=-\mathrm{0.0040455435},\\ C_{\mathrm{21}}=-\mathrm{0.0602863302},\\ C_{\mathrm{22}}=\mathrm{0.7021897011},\\ C_{\mathrm{23}}=-\mathrm{1.1196100026}.\end{array}$$Substituting the values of $C_{11}$, $C_{12}$, $C_{13}$, and $C_{21}$, $C_{22}$, and $C_{23}$, into (37) and (38), one obtains$$\begin{array}{c}\text{(40)}& u_{\mathrm{1}}\left(t\right)=-\left(\mathrm{0.1526946}{e}^t-\mathrm{0.04241517}\right)t^{\mathrm{6}}+\left(\mathrm{1.488772}{e}^t-\mathrm{0.2163173}\right)t^{\mathrm{5}}-\left(\mathrm{9.47979}{e}^t+\mathrm{0.1627504}\right)t^{\mathrm{4}}+\left(\mathrm{42.16258}{e}^t-\mathrm{1.526946}{e}^{-t}+\mathrm{1.162051}\right)t^{\mathrm{3}}+\left(-\mathrm{132.3468}{e}^t-\mathrm{7.63473}{e}^{-t}+\mathrm{3.832068}\right)t^{\mathrm{2}}+\left(\mathrm{277.0864}{e}^t-\mathrm{16.79641}{e}^{-t}+\mathrm{22.74883}\right)t+\left(\mathrm{1.017005}\left(t-\mathrm{1}\right)-\mathrm{298.8182}+\mathrm{0.5085023}\left(t^{\mathrm{3}}-\mathrm{3}{t}^{\mathrm{2}}+\mathrm{6}t-\mathrm{6}\right)\right)e^t-\mathrm{16.79641}{e}^{-t}+\mathrm{320.6826},u_{\mathrm{2}}\left(t\right)=\mathrm{0.000003651786}{t}^{\mathrm{6}}+\left(-\mathrm{0.000009129464}{e}^{-t}+\mathrm{0.003148389}\right)t^{\mathrm{5}}+\left(\mathrm{0.005375128}{e}^{-t}-\mathrm{0.01171288}\right)t^{\mathrm{4}}+\left(\mathrm{0.0216831}{e}^{-t}-\mathrm{0.1448423}\right)t^{\mathrm{3}}+\left(\mathrm{0.00715933}{e}^{-\mathrm{2}t}+\mathrm{0.03383929}{e}^{-t}+\mathrm{0.43758}\right)t^{\mathrm{2}}+\left(\mathrm{0.01442821}{e}^{-\mathrm{2}t}-\mathrm{0.1490952}{e}^{-t}-\mathrm{1.197225}\right)t+\left(-\mathrm{1.179469}\left(t+\mathrm{1}\right)+\mathrm{0.0003527703}\left(e^t\left(\mathrm{12}{t}^{\mathrm{5}}-\mathrm{45}{t}^{\mathrm{4}}+\mathrm{160}{t}^{\mathrm{3}}-\mathrm{240}{t}^{\mathrm{3}}-\mathrm{240}t+\mathrm{720}\right)+\mathrm{20}{t}^{\mathrm{4}}+\mathrm{80}{t}^{\mathrm{3}}+\mathrm{120}{t}^{\mathrm{2}}-\mathrm{480}t-\mathrm{720}\right)-\mathrm{0.559805}\left(t^{\mathrm{3}}+\mathrm{3}{t}^{\mathrm{2}}+\mathrm{6}t+\mathrm{6}\right)-\mathrm{0.2269294}\right)e^{-t}+\mathrm{0.00732366}{e}^{-\mathrm{2}t}+\mathrm{5.757905}.\end{array}$$To derive a solution to (34), for $0\le t\le 1$, by MOHAM, we consider the following initial approximation:$$\begin{array}{}\text{(41)}& u_{i\mathrm{0}}^j\left(t_j\right)={\alpha }_i^j,\hspace{1em}i=\mathrm{1,2},\hspace{0.17em}\hspace{0.17em}j=\mathrm{1},\dots ,\mathrm{5}.\end{array}$$Now, consider the auxiliary function $H_i\left(P,t\right)$, as the following:$$\begin{array}{}\text{(42)}& H_i\left(P,t\right)=\left(C_{i\mathrm{1}}^j+C_{i\mathrm{2}}^{j}t+C_{i\mathrm{3}}^j{t}^{\mathrm{2}}\right)P.\end{array}$$where $C_{iq}^j$, $i=\mathrm{1,2}$, $j=\mathrm{1},\dots ,\mathrm{5}$, $q=\mathrm{1,2},\mathrm{3}$, are still unknown.

Regarding (17), the first-order approximate solution, for $m=1$, is as follows: $$\begin{array}{}\text{(43)}& {\tilde{u}}_i^j\left(t\right)=u_{i\mathrm{0}}^j\left(t\right)+u_{i\mathrm{1}}^j\left(t\right),\hspace{1em}i=\mathrm{1,2},\hspace{0.17em}\hspace{0.17em}j=\mathrm{1},\dots ,\mathrm{5},\end{array}$$where$$\begin{array}{}\text{(44)}& u_{\mathrm{11}}^j\left(t\right)=\left(C_{\mathrm{11}}^j+C_{\mathrm{12}}^{j}t+C_{\mathrm{13}}^j{t}^{\mathrm{2}}\right)\left(u_{\mathrm{10}}^j\left(t\right)-\mathrm{1}+\frac{t^{\mathrm{2}}}{\mathrm{2}}-{{\displaystyle \int }}_{\mathrm{0}}^t\left(u_{\mathrm{10}}^j\left(s\right)+s{e}^s{u}_{\mathrm{20}}^j\left(s\right)\right)ds\right),\text{(45)}& u_{\mathrm{21}}^j\left(t\right)=\left(C_{\mathrm{21}}^j+C_{\mathrm{22}}^{j}t+C_{\mathrm{23}}^j{t}^{\mathrm{2}}\right)\left(u_{\mathrm{20}}^j\left(t\right)-\mathrm{1}-\frac{t^{\mathrm{2}}}{\mathrm{2}}+{{\displaystyle \int }}_{\mathrm{0}}^t\left(s{e}^s{u}_{\mathrm{10}}^j\left(s\right)+u_{\mathrm{20}}^j\left(s\right)\right)ds\right),\end{array}$$Substitution of (44) and (45) into (43), regarding (41) and using (18)–(20), determines the residual and the functional $J_1^j\left(C_{1q}^j\right)$, and $J_2^j\left(C_{2q}^j\right)$ for $q=\mathrm{1,2},3$, respectively$$\begin{array}{c}\text{(46)}& R_{\mathrm{1}}^j\left(t,C_{\mathrm{1}q}^j\right)={\tilde{u}}_{\mathrm{1}}^j\left(t\right)-\mathrm{1}+\frac{t^{\mathrm{2}}}{\mathrm{2}}-{{\displaystyle \int }}_{\mathrm{0}}^t\left(u_{\mathrm{1}}^j\left(s\right)+s{e}^s{u}_{\mathrm{2}}^j\left(s\right)\right)ds,R_{\mathrm{2}}^j\left(t,C_{\mathrm{2}q}^j\right)={\tilde{u}}_{\mathrm{2}}^j\left(t\right)-\mathrm{1}-\frac{t^{\mathrm{2}}}{\mathrm{2}}+{{\displaystyle \int }}_{\mathrm{0}}^t\left(s{e}^s{u}_{\mathrm{1}}^j\left(s\right)+u_{\mathrm{2}}^j\left(s\right)\right)ds,\\ J_{\mathrm{1}}^j\left(C_{\mathrm{1}q}^j\right)={{\displaystyle \int }}_{t_{\gamma }}^{t_{\gamma }+h}{\left(R_{\mathrm{1}}^j\right)}^{\mathrm{2}}\left(s,C_{\mathrm{1}q}^j\right)ds,\\ J_{\mathrm{2}}^j\left(C_{\mathrm{2}q}^j\right)={{\displaystyle \int }}_{t_{\gamma }}^{t_{\gamma }+h}{\left(R_{\mathrm{2}}^j\right)}^{\mathrm{2}}\left(s,C_{\mathrm{2}q}^j\right)ds.\end{array}$$Minimization condition required$$\begin{array}{c}\text{(47)}& \frac{\partial J_{\mathrm{1}}^j}{\partial C_{\mathrm{11}}^j}=\frac{\partial J_{\mathrm{1}}^j}{\partial C_{\mathrm{12}}^j}=\frac{\partial J_{\mathrm{1}}^j}{\partial C_{\mathrm{13}}^j}=\mathrm{0},\frac{\partial J_{\mathrm{2}}^j}{\partial C_{\mathrm{21}}^j}=\frac{\partial J_{\mathrm{2}}^j}{\partial C_{\mathrm{22}}^j}=\frac{\partial J_{\mathrm{2}}^j}{\partial C_{\mathrm{23}}^j}=\mathrm{0}.\end{array}$$Parameters that control the convergence, $C_{iq}^j$, $i=\mathrm{1,2}$, $j=\mathrm{1},\dots ,\mathrm{5}$, $q=\mathrm{1,2},\mathrm{3}$, are presented in Table 1, by considering ${\alpha }_{\mathrm{1}}^{\mathrm{1}}=\mathrm{0}$, ${\alpha }_{\mathrm{2}}^{\mathrm{1}}=\mathrm{0}$, and $h=\mathrm{0.2}$, for $t_{\mathrm{0}}=\mathrm{0}$ up to $t_{\mathrm{5}}=T=\mathrm{1}.$

By substituting the values of the control parameters $C_{iq}^j$’s into (44) and (45), one obtains$$\begin{array}{c}\text{(48)}& {\tilde{u}}_{\mathrm{1}}\left(t\right)=\left\{\begin{array}{cc}-\left(\mathrm{0.5}{t}^{\mathrm{2}}-\mathrm{1}\right)\left(\mathrm{1.5241}{t}^{\mathrm{2}}+\mathrm{0.95441}t+\mathrm{1.0008}\right),\hfill & \mathrm{0}\le t<\mathrm{0.2},\hfill \\ \left(\mathrm{0.083435}{t}^{\mathrm{2}}+\mathrm{0.28081}t+\mathrm{0.95613}\right)\left(\mathrm{0.14906}{t}^{\mathrm{5}}{e}^t+\mathrm{1.259}{t}^{\mathrm{4}}{e}^t+\mathrm{5.8342}{t}^{\mathrm{3}}{e}^t\right.\hfill & \hfill \\ -\mathrm{18.53}{t}^{\mathrm{2}}{e}^t-\mathrm{38.061}{e}^t+\mathrm{1.0008}t+\left(\mathrm{0.5}{t}^{\mathrm{2}}-\mathrm{1}\right)\left(\mathrm{1.5241}{t}^{\mathrm{2}}+\mathrm{0.95441}t+\mathrm{1.0008}\right)\hfill & \hfill \\ \left.-\mathrm{0.15241}{t}^{\mathrm{5}}-\mathrm{0.1193}{t}^{\mathrm{4}}+\mathrm{0.34122}{t}^{\mathrm{3}}-\mathrm{0.022794}{t}^{\mathrm{2}}+\mathrm{38.061}t{e}^t+\mathrm{39.061}\right)\hfill & \hfill \\ -\left(\mathrm{0.5}{t}^{\mathrm{2}}-\mathrm{1}\right)\left(\mathrm{1.5}{\mathrm{241}t}^{\mathrm{2}}+\mathrm{0.95441}t+\mathrm{1.0008}\right),\hfill & \mathrm{0.2}\le t<\mathrm{0.4},\hfill \\ ⋮\hfill & \hfill \end{array}\right.{\tilde{u}}_{\mathrm{2}}\left(t\right)=\left\{\begin{array}{cc}\left(\mathrm{0.5}{t}^{\mathrm{2}}+\mathrm{1}\right)\left(\mathrm{0.29812}{t}^{\mathrm{2}}-\mathrm{1.0273}t+\mathrm{1.0006}\right),\hfill & \mathrm{0}\le t<\mathrm{0.2},\hfill \\ -\left(-\mathrm{0.39731}{t}^{\mathrm{2}}+\mathrm{0.061232}t+\mathrm{0.91564}\right)\left(\left(\mathrm{0.5}{t}^{\mathrm{2}}+\mathrm{1}\right)\left(\mathrm{0.29812}{t}^{\mathrm{2}}-\mathrm{1.0273}t\right.\right.\hfill & \hfill \\ \left.+\mathrm{1.0006}\right)-\mathrm{0.5}{t}^{\mathrm{2}}+\mathrm{1.8504}{e}^{-\left(t+\mathrm{18}\right)}\left(e^{\left(t+\mathrm{16}\right)}\left(\mathrm{1.6111}{t}^{\mathrm{5}}-\mathrm{6.94}{t}^{\mathrm{4}}\right)\right.\hfill & \hfill \\ -\mathrm{5.0717}{e}^{\left(t+\mathrm{19}\right)}+\mathrm{4.1183}{t}^{\mathrm{5}}{e}^{\mathrm{17}}+e^{\mathrm{18}}\left(\mathrm{2.317}{t}^{\mathrm{4}}+\mathrm{8.715}{t}^{\mathrm{3}}\right)\hfill & \hfill \\ \left.\left.+e^{\left(t+\mathrm{17}\right)}\left(\mathrm{1.4383}{t}^{\mathrm{3}}-\mathrm{2.776}{t}^{\mathrm{2}}+\mathrm{5.4076}t\right)+e^{\mathrm{19}}\left(\mathrm{2.5629}{t}^{\mathrm{2}}+\mathrm{5.0717}t+\mathrm{5.0717}\right)\right)-1\right)\hfill & \hfill \\ +\left(\mathrm{0.5}{t}^{\mathrm{2}}+\mathrm{1}\right)\left(\mathrm{0.29812}{t}^{\mathrm{2}}-\mathrm{1.0273}t+\mathrm{1.0006}\right),\hfill & \mathrm{0.2}\le t<\mathrm{0.4},\hfill \\ ⋮\hfill & \hfill \end{array}\right.\end{array}$$Expressions of the solution of the other interval are too long, that is why they are not presented here. Interested readers can calculate the solutions by the program appeared in Appendix A.

First-order MOHAM approximate solutions and three-order OHAM approximate solutions can be compared with exact solutions in Tables 2 and 3, and plots are presented in Figures 1 and 2.

Absolute errors for OHAM and MOHAM are plotted in Figures 3 and 4.

Exact $u_i$, OHAM $u_i$, and MOHAM $u_i$ stand for exact solution of $u_i$, the solution of $u_i$ by OHAM, and solution of by MOHAM, respectively.

Table 1

Values of the control parameters $C_{ij}$.

$$\mathrm{j}$$	$$C_{11}^j$$	$$C_{12}^j$$	$$C_{13}^j$$	$$C_{21}^j$$	$$C_{22}^j$$	$$C_{23}^j$$
1	-1.00082	-0.95441	-1.52407	-1.00061	1.02732	-0.29811
2	-0.95613	-0.28080	-0.08343	-0.91564	-0.06123	0.39730
3	-0.77599	-0.63069	-0.00753	0.00000	0.00000	0.00000
4	0.00000	0.00000	0.00000	-0.36062	-0.91215	0.19472
5	0.00000	0.00000	0.00000	-0.89126	-0.01338	0.00040

Table 2

The results of applying OHAM, MOHAM, and the exact solution, for Example 1.

$${\mathrm{t}}_{\mathrm{j}}$$	Exact$u_1$	Exact$u_2$	OHAM$u_1$	OHAM$u_2$	MOHAM$u_1$	MOHAM$u_2$
0.0	1.000000	1.000000	1.000000	1.000000	1.000822	1.000615
0.1	1.105170	0.904837	1.106625	0.898382	1.105964	0.905368
0.2	1.221402	0.818730	1.225525	0.815433	1.227614	0.823217
0.3	1.349858	0.740818	1.360566	0.748515	1.350942	0.740830
0.4	1.491824	0.670320	1.515929	0.695399	1.494814	0.670041
0.5	1.648721	0.606530	1.696140	0.654166	1.648851	0.606287
0.6	1.822118	0.548811	1.906111	0.623127	1.821875	0.549821
0.7	2.013752	0.496585	2.151164	0.600764	2.012654	0.496810
0.8	2.225540	0.449328	2.437011	0.585686	2.223236	0.448796
0.9	2.459603	0.406569	2.769628	0.576594	2.456573	0.406549
1.0	2.718281	0.367879	3.154950	0.572265	2.717034	0.367906

Table 3

Absolute errors of applying OHAM and MOHAM for Example 1.

$$t_j$$	Ab. Error(OHAM$u_1$)	Ab. Error(OHAM$u_2$)	Ab. Error(MOHAM$u_1$)	Ab. Error(MOHAM$u_2$)
0.0	0.000000	0.000000	0.000822	0.000615
0.1	0.001455	0.006455	0.000775	0.000531
0.2	0.004123	0.003296	0.006211	0.004487
0.3	0.010708	0.007697	0.001083	0.004486
0.4	0.024105	0.025079	0.002990	1.2469${\mathrm{e}}^{-\mathrm{5}}$
0.5	0.047419	0.047635	0.000130	0.000242
0.6	0.083993	0.074315	0.000242	0.001010
0.7	0.137412	0.104179	0.001098	0.000225
0.8	0.211471	0.136357	0.002304	0.000532
0.9	0.310025	0.170024	0.003030	2.0620${\mathrm{e}}^{-\mathrm{5}}$
1.0	0.436669	0.204386	0.001247	2.7127${\mathrm{e}}^{-\mathrm{5}}$

[figure omitted; refer to PDF] [figure omitted; refer to PDF] [figure omitted; refer to PDF] [figure omitted; refer to PDF]

Example 2.

let us consider the following nonlinear system of Volterra integral equations of the second kind, with the exact solutions $u_1\left(t\right)=t$ and $u_2\left(t\right)=t$. $$\begin{array}{c}\text{(49)}& u_{\mathrm{1}}\left(t\right)=t-t^{\mathrm{2}}-{{\displaystyle \int }}_{\mathrm{0}}^t\left(u_{\mathrm{1}}\left(s\right)+u_{\mathrm{2}}\left(s\right)\right)ds,u_{\mathrm{2}}\left(t\right)=t-\frac{t^{\mathrm{2}}}{\mathrm{2}}-\frac{t^{\mathrm{3}}}{\mathrm{3}}-{{\displaystyle \int }}_{\mathrm{0}}^t\left({u_{\mathrm{1}}}^{\mathrm{2}}\left(s\right)+u_{\mathrm{2}}\left(s\right)\right)ds.\end{array}$$Following the procedure in Example 1, we have$$\begin{array}{c}\text{(50)}& C_{\mathrm{11}}=\mathrm{0},C_{\mathrm{12}}=-\mathrm{0.0000000008},\\ C_{\mathrm{13}}=-\mathrm{1.9560359039},\\ C_{\mathrm{21}}=\mathrm{0.0000000018},\\ C_{\mathrm{22}}=-\mathrm{0.0000000033},\\ C_{\mathrm{23}}=-\mathrm{1.4774464040},\\ u_{\mathrm{1}}\left(t\right)=\mathrm{8.436182}{e}^{-\mathrm{35}}{t}^{\mathrm{8}}+\mathrm{8.3237}{e}^{-\mathrm{34}}{t}^{\mathrm{7}}+\mathrm{8.756976}{e}^{-\mathrm{22}}{t}^{\mathrm{6}}-e^{-\mathrm{21}}\left(\mathrm{3.50279}{t}^{\mathrm{5}}+\mathrm{2.189244}{t}^{\mathrm{4}}\right)-\mathrm{4.533431}{e}^{-\mathrm{22}}{t}^{\mathrm{3}}\left(-\mathrm{2}{t}^{\mathrm{3}}+\mathrm{8}{t}^{\mathrm{2}}+\mathrm{5}t-\mathrm{30}\right)+\mathrm{1.313546}{e}^{-\mathrm{20}}{t}^{\mathrm{3}}-\mathrm{0.163003}{t}^{\mathrm{2}}\left(t^{\mathrm{2}}+\mathrm{6}t-\mathrm{12}\right)-t^{\mathrm{2}}+t,\\ u_{\mathrm{2}}\left(t\right)=\mathrm{3.777777}{e}^{-\mathrm{29}}{t}^{\mathrm{9}}-\mathrm{3.134551}{e}^{-\mathrm{28}}{t}^{\mathrm{8}}-\mathrm{6.437708}{e}^{-\mathrm{20}}{t}^{\mathrm{7}}-\mathrm{4.656609}{e}^{-\mathrm{19}}{t}^{\mathrm{6}}+e^{-\mathrm{18}}\left(\mathrm{1.93775}{t}^{\mathrm{5}}-\mathrm{1.464579}{t}^{\mathrm{4}}\right)+\mathrm{1.922551}{e}^{-\mathrm{22}}{t}^{\mathrm{3}}\left(\mathrm{60}{t}^{\mathrm{4}}+\mathrm{434}{t}^{\mathrm{3}}-\mathrm{1806}{t}^{\mathrm{2}}+\mathrm{1365}t+\mathrm{420}\right)+\mathrm{0.02462411}{t}^{\mathrm{2}}\left(\mathrm{12}{t}^{\mathrm{3}}-\mathrm{35}{t}^{\mathrm{2}}+\mathrm{10}t+\mathrm{30}\right)-\mathrm{0.3333333}{t}^{\mathrm{3}}-\mathrm{0.5}{t}^{\mathrm{2}}+t,\end{array}$$Parameters that control the convergence, $C_{iq}^j$, $i=\mathrm{1,2}$, $j=\mathrm{1},\dots ,\mathrm{5}$, $q=\mathrm{1,2},\mathrm{3}$ are presented in Table 4, by considering ${\alpha }_{\mathrm{1}}^{\mathrm{1}}=\mathrm{0}$, ${\alpha }_{\mathrm{2}}^{\mathrm{1}}=\mathrm{0}$, $h=\mathrm{0.2}$, for $t_{\mathrm{0}}=\mathrm{0}$ up to $t_{\mathrm{5}}=T=\mathrm{1}.$

Therefore, the approximate solutions can be written in the following form:$$\begin{array}{c}\text{(51)}& {\tilde{u}}_{\mathrm{1}}\left(t\right)=\left\{\begin{array}{cc}\left(t-t^{\mathrm{2}}\right)\left(\mathrm{1.83126}{t}^{\mathrm{2}}+\mathrm{0.916398}t\mathrm{1.00294}\right),\hfill & \mathrm{0}\le t<\mathrm{0.2},\hfill \\ \left(t-t^{\mathrm{2}}\right)\left(\mathrm{1.83126}{t}^{\mathrm{2}}+\mathrm{0.916398}t\mathrm{1.00294}\right),\hfill & \mathrm{0.2}\le t<\mathrm{0.4},\hfill \\ -\left(\mathrm{0.335879}{t}^{\mathrm{2}}+\mathrm{0.39467}t+\mathrm{0.855231}\right)\left(\mathrm{1.54198}{e}^{-\mathrm{19}}{t}^{\mathrm{2}}\left(e^{\mathrm{18}}\left({\mathrm{3.03659}t}^{\mathrm{3}}\right.\right.\right.\hfill & \hfill \\ \left.\left.-\mathrm{1.67635}{t}^{\mathrm{2}}-\mathrm{6.49617}\right)+e^{\mathrm{17}}\left(\mathrm{2.50648}{t}^{\mathrm{4}}+\mathrm{2.17358}t\right)\right)+\left(t-t^{\mathrm{2}}\right)\hfill & \hfill \\ .\left.\left(\mathrm{1.83126}{t}^{\mathrm{2}}+\mathrm{0.916398}t+\mathrm{1.00294}\right)+t^{\mathrm{2}}-t\right)+\left(t-t^{\mathrm{2}}\right)\left(\mathrm{1.83126}{t}^{\mathrm{2}}\right.\hfill & \hfill \\ \left.+\mathrm{0.916398}t+\mathrm{1.00294}\right),\hfill & \mathrm{0.4}\le t<\mathrm{0.6},\hfill \\ ⋮\hfill & \hfill \end{array}\right.{\tilde{u}}_{\mathrm{2}}\left(t\right)=\left\{\begin{array}{cc}-\left(\mathrm{0.486219}t+\mathrm{0.695688}{t}^{\mathrm{2}}+\mathrm{1.00045}\right)\left(\mathrm{0.5}{t}^{\mathrm{2}}-\mathrm{t}+\mathrm{0.333333}{t}^{\mathrm{3}}\right)\mathbf{,}\hfill & \mathrm{0}\le t<\mathrm{0.2},\hfill \\ -\left(\mathrm{0.486219}t+\mathrm{0.695688}{t}^{\mathrm{2}}+\mathrm{1.00045}\right)\left(\mathrm{0.5}{t}^{\mathrm{2}}-\mathrm{t}+\mathrm{0.333333}{t}^{\mathrm{3}}\right)\mathbf{,}\hfill & \mathrm{0.2}\le t<\mathrm{0.4},\hfill \\ \left(\mathrm{1.88661}{t}^{\mathrm{2}}-\mathrm{0.994006}t+\mathrm{1.14645}\right)\left(t+\mathrm{2.7391}{e}^{-\mathrm{33}}{t}^{\mathrm{2}}\left(\mathrm{1.36035}{e}^{\mathrm{32}}{t}^{\mathrm{7}}-\mathrm{1.52911}{e}^{\mathrm{32}}{t}^{\mathrm{6}}\right.\right.\hfill & \hfill \\ +\mathrm{6.01834}{e}^{\mathrm{31}}{t}^{\mathrm{5}}-\mathrm{2.47255}{e}^{\mathrm{32}}{t}^{\mathrm{4}}+\mathrm{9.73078}{e}^{\mathrm{31}}{t}^{\mathrm{3}}-\mathrm{4.97417}{e}^{\mathrm{30}}{t}^{\mathrm{2}}+\mathrm{1.20707}{e}^{\mathrm{32}}t\hfill & \hfill \\ \left.+\mathrm{1.82623}{e}^{\mathrm{32}}\right)-\mathrm{0.333333}{t}^{\mathrm{3}}-\mathrm{0.5}{t}^{\mathrm{2}}+\left(\mathrm{0.695688}{t}^{\mathrm{2}}+\mathrm{0.486219}t+\mathrm{1.00045}\right)\hfill & \hfill \\ .\left.\left(\mathrm{0.333333}{t}^{\mathrm{3}}+\mathrm{0.5}{t}^{\mathrm{2}}-t\right)\right)-\left(\mathrm{0.695688}{t}^{\mathrm{2}}+\mathrm{0.486219}t+\mathrm{1.00045}\right)\hfill & \hfill \\ .\left(\mathrm{0.333333}{t}^{\mathrm{3}}+\mathrm{0.5}{t}^{\mathrm{2}}-t\right),\hfill & 0.4\le t<0.6,\hfill \\ ⋮\hfill & \hfill \end{array}\right.\end{array}$$Due to the long expressions of the solutions on other intervals, the results are not reported, but related programs will be appeared in Appendix B.

First-order MOHAM approximate solution and three-order OHAM approximate solution can be compared with an exact solution in Tables 5 and 6, and plots are presented in Figures 5 and 6.

Absolute errors for OHAM and MOHAM are plotted in Figures 7 and 8.

Table 4

Values of the control parameters $C_{ij}$.

$$\mathrm{j}$$	$$C_{11}^j$$	$$C_{12}^j$$	$$C_{13}^j$$	$$C_{21}^j$$	$$C_{22}^j$$	$$C_{23}^j$$
1	-1.002941	-0.916397	-1.831260	-1.000448	-0.486218	-0.695688
2	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000
3	-0.855231	-0.394669	-0.335878	-1.146447	0.994005	-1.886613
4	-0.771655	0.098153	-0.628421	0.000000	0.000000	0.000000
5	-0.085517	-0.993931	0.069185	-0.666758	0.328380	-0.669028

Table 5

The results of applying OHAM, MOHAM, and the exact solution for Example 2.

$${\mathrm{t}}_{\mathrm{j}}$$	Exact$u_{\mathrm{1}}$	Exact$u_{\mathrm{2}}$	OHAM$u_{\mathrm{1}}$	OHAM$u_{\mathrm{2}}$	MOHAM$u_{\mathrm{1}}$	MOHAM$u_{\mathrm{2}}$
0.00.10.20.30.40.50.60.70.80.91.0	0.00.10.20.30.40.50.60.70.80.91.0	0.00.10.20.30.40.50.60.70.80.91.0	0.0000000.1085660.2301560.3583160.4861990.6065690.7117950.7938600.8443510.8544670.815014	0.0000000.1022160.2075670.3128700.4135840.5041630.5784100.6298350.6520060.6389030.585276	0.0000000.1001600.2015150.3029610.3990000.4995310.6009590.7002960.7996660.8997601.000794	0.0000000.0999700.1995920.2973950.3901320.4999590.5999660.7002290.8009620.9004600.999225

Table 6

Absolute errors of applying OHAM and MOHAM for Example 2.

$t_j$	Ab. Error(OHAM$u_{\mathrm{1}}$)	Ab. Error(OHAM$u_{\mathrm{2}}$)	Ab. Error(MOHAM$u_{\mathrm{1}}$)	Ab. Error(MOHAM$u_{\mathrm{2}}$)
0.0	0.000000	0.000000	0.000000	0.000000
0.1	0.008566	0.002216	0.000160	3.0000$e^{-5}$
0.2	0.030156	0.007567	0.001515	0.000408
0.3	0.058316	0.012870	0.002961	0.002605
0.4	0.086199	0.013584	0.001000	0.009868
0.5	0.106569	0.004163	0.000469	4.1000$e^{-5}$
0.6	0.111795	0.021589	0.000959	3.4000$e^{-5}$
0.7	0.093860	0.070165	0.000296	0.000229
0.8	0.044351	0.147994	0.000334	0.000962
0.9	0.045533	0.261097	0.000240	0.000460
1.0	0.184986	0.414724	0.000794	0.000775

[figure omitted; refer to PDF] [figure omitted; refer to PDF] [figure omitted; refer to PDF] [figure omitted; refer to PDF]

7. Conclusion and Discussion

In this paper, two well known approaches, OHAM, and MOHAM, have been applied for solving linear and nonlinear systems of Volterra integral equations of the second kind. The results of applying these two approaches are presented in Tables and are plotted in Figures. Tables 3 and 6 show the absolute errors of applying OHAM and MOHAM at some selected points. Comparing numerical results, reveal that MOHAM is more accurate than OHAM, especially for the points farther from the initial point. Moreover MOHAM is very efficient and convenient to use for finding approximation solution for system of Volterra integral equations of the second kind.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of the paper.

Appendix

Appendix

A. Example 1

 

clear; clc; format long g; close all

 

T=1; dt=0.02; tt=0:dt:T; tic

 

%$\ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast$

 

h=0.2;

 

%$\ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast$

 

u1exact=@(x)exp(x);

 

u2exact=@(x)exp(-x);

 

ttt=0:h:T; alfa1=0; alfa2=0; v=5;

 

options = optimoptions('fsolve', 'Display', 'off',

 

'Algorithm', 'trust-region-dogleg',...

 

'MaxFunEvals',1e6,'MaxIter',1000,'TolFun',1e-10,'TolX',1e-10);

 

syms x t c1 c2 c3

 

u10=alfa1; u20=alfa2;

 

H1=c1+ c2$\ast$x+c3$\ast {\mathrm{x}}^{\wedge }$2;

 

h1=1-(${\mathrm{x}}^{\wedge }$2)/2; h2=1+(${\mathrm{x}}^{\wedge }$2)/2;

 

for i=1:length(ttt)-1

 

f10=subs(u10,x,t); f20=t$\ast$exp(t)$\ast$subs(u20,x,t); f30=t$\ast$exp(-t)$\ast$subs(u10,x,t); f40=subs(u20,x,t);

 

f1=f10+f20;

 

f2=f30+f40;

 

in1=int(f1,t,0,x );

 

u11=H1$\ast$ ( u10 -h1 -in1 );

 

in2=int(f2,t,0,x );

 

u21=H1$\ast$ ( u20 -h2 +in2 );

 

ut1=u10+u11; ut2=u20+u21;

 

g1=subs(ut1,x,t); g2=subs(ut2,x,t);

 

ig1=g1+t$\ast$exp(t)$\ast$g2; int1=int(ig1,t,0,x);

 

R1=ut1-h1-int1;

 

I1=vpa(R$1^{\wedge }$2,v);

 

I11=int(I1,x,ttt(i),ttt(i+1));

 

J1=vpa( I11, v);

 

fprintf('∖n∖n J1(x) = %s, ∖n',char((vpa(J1,v))) );

 

ig2=t$\ast$exp-t)$\ast$g1+g2; int2=int(ig2,t,0,x);

 

R2=ut2-h2+int2;

 

I2=vpa(R$2^{\wedge }$2,v);

 

I21=int(I2,x,ttt(i),ttt(i+1));

 

J2=vpa( I21, v);

 

fprintf('∖n∖n J2(x) = %s, ∖n',char((vpa(J2,v))) );

 

J11 = vpa(diff(J1,c1),v);

 

J12 = vpa(diff(J1,c2),v);

 

J13 = vpa(diff(J1,c3),v);

 

J21 = vpa(diff(J2,c1),v);

 

J22 = vpa(diff(J2,c2),v);

 

J23 = vpa(diff(J2,c3),v);

 

global dJ11 dJ12 dJ13 dJ21 dJ22 dJ23;

 

dJ11=@(c1,c2,c3)eval(char(J11)); dJ12=@(c1,c2,c3)eval(char(J12)); dJ13=@(c1,c2,c3)eval(char(J13));

 

dJ21=@(c1,c2,c3)eval(char(J21)); dJ22=@(c1,c2,c3)eval(char(J22)); dJ23=@(c1,c2,c3)eval(char(J23));

 

options = optimoptions('fsolve', 'Display', 'off',

 

'Algorithm', 'trust-region-dogleg',...

 

'MaxFunEvals',1e6,'MaxIter',1000,'TolFun',1e-10,'TolX',1e-10);

 

cc1=fsolve(@Sifun1,$[0\hspace{0.17em}0\hspace{0.17em}0]$,options);

 

cc2=fsolve(@Sifun2,$[0\hspace{0.17em}0\hspace{0.17em}0]$,options);

 

fprintf('∖n∖n∖n∖t∖t∖t∖t∖t∖tc11 = %1.10f,∖n∖n ∖t∖t∖t∖t∖t∖tc12 = %1.10f,∖n∖n ∖t∖t∖t∖t∖t∖tc13 = %1.10f,∖n∖n',cc1(1),cc1(2),cc1(3));

 

fprintf('∖n∖n∖n∖t∖t∖t∖t∖t∖tc21 = %1.10f,∖n∖n ∖t∖t∖t∖t∖t∖tc22 = %1.10f,∖n∖n ∖t∖t∖t∖t∖t∖tc23 = %1.10f, ∖n∖n',cc2(1),cc2(2),cc2(3));

 

u1=( subs( ut1,[c1,c2,c3],[cc1(1),cc1(2),cc1(3)] ) );

 

u2=( subs( ut2,[c1,c2,c3],[cc2(1),cc2(2),cc2(3)] ) );

 

fprintf('∖n∖n u1(%g) = %s, ∖n',ttt(i),char(vpa(u1,v)) );

 

fprintf('∖n∖n u2(%g) = %s, ∖n',ttt(i),char(vpa(u2,v)) );

 

yy1=@(x)eval(char(u1));

 

yy2=@(x)eval(char(u2));

 

yy=@(x)eval(char(u));

 

mt=ttt(i):dt:ttt(i+1);

 

for j=1:length(mt)

 

yyn1(j)=yy1(mt(j));

 

yyn2(j)=yy2(mt(j));

 

y1exn(j)=u1exact(mt(j));

 

y2exn(j)=u2exact(mt(j));

 

end

 

yapp$\{\mathrm{1},\mathrm{i}\}$=yyn1; yapp$\{\mathrm{2},\mathrm{i}\}$=yyn2; yex$\{\mathrm{1},\mathrm{i}\}$=y1exn; yex$\{\mathrm{2},\mathrm{i}\}$=y2exn;

 

u10=u1; u20=u2;

 

end

 

yapprox1=yapp$\{\mathrm{1,1}\}$; yapprox2=yapp$\{\mathrm{2,1}\}$; yexact1=yex$\{\mathrm{1,1}\}$; yexact2=yex$\{\mathrm{2,1}\}$;

 

for i=2:length(ttt)-1

 

yapprox1=[yapprox1,yapp$\{\mathrm{1},\mathrm{i}\}$(2:end)];

 

yapprox2=[yapprox2,yapp$\{\mathrm{2},\mathrm{i}\}$(2:end)];

 

yexact1=[yexact1,yex$\{\mathrm{1},\mathrm{i}\}$(2:end)];

 

yexact2=[yexact2,yex$\{\mathrm{2},\mathrm{i}\}$(2:end)];

 

end

 

th=0:0.1:T; er1=abs(yexact1-yapprox1); er2=abs(yexact2-yapprox2);

 

for i=1:length(th)

 

k(i)=round(th(i)/dt)+1;

 

end

 

L=$\left\{\hspace{0.17em}\right.$'t', 'Exact1', 'Exact2', ' MOHAM1', ' MOHAM2', 'Error1', 'Error2'$\left.\hspace{0.17em}\right\}$;

 

table(th',yexact1(k)', yexact2(k)', yapprox1(k)', yapprox2(k)',er1(k)',er2(k) ','VariableNames',L)

 

figure('Units','characters','Name','p1','Position',

 

$[25\hspace{0.17em}\hspace{0.17em}6\hspace{0.17em}\hspace{0.17em}160\hspace{0.17em}\hspace{0.17em}45]$ );

 

plot(tt,yexact1,'linewidth',2);hold on;

 

plot(tt,yapprox1,'o','MarkerSize',8); grid on; grid minor; axis tight;

 

xlabel('x','fontweight','Bold');

 

legend('u_1 Exact solution','u_1 MOHAM solution', 'Location','northwest');

 

set( gca, 'fontname', 'Euclid', 'fontsize',12 );

 

figure('Units','characters','Name', 'p2', 'Position',

 

$[25\hspace{0.17em}\hspace{0.17em}6\hspace{0.17em}\hspace{0.17em}160\hspace{0.17em}\hspace{0.17em}45]$ );

 

plot(tt,yexact2,'linewidth',2);hold on;

 

plot(tt,yapprox2,'o','MarkerSize',8); grid on; grid minor; axis tight;

 

xlabel('x','fontweight','Bold');

 

legend('u_2 Exact solution','u_2 MOHAM solution', 'Location', 'northwest');

 

set( gca, 'fontname', 'Euclid', 'fontsize',12 );

 

figure('Units','characters','Name','Error1', 'Position',

 

$[25\hspace{0.17em}\hspace{0.17em}6\hspace{0.17em}\hspace{0.17em}160\hspace{0.17em}\hspace{0.17em}45]$ );

 

plot(tt,er1,'linewidth',2);

 

grid on; grid minor; axis tight; xlabel('x','fontweight', 'Bold');

 

ch=ylabel('u_1 Absolute Error', 'fontweight', 'Bold');

 

set( gca, 'fontname', 'Euclid','fontsize',12 );

 

figure('Units','characters','Name','Error2', 'Position',$[25\hspace{0.17em}\hspace{0.17em}6\hspace{0.17em}\hspace{0.17em}160\hspace{0.17em}\hspace{0.17em}45]$);

 

plot(tt,er2,'linewidth',2);

 

grid on; grid minor; axis tight; xlabel('x', 'fontweight', 'Bold');

 

ch=ylabel('u_1 Absolute Error','fontweight','Bold');

 

set( gca, 'fontname', 'Euclid', 'fontsize',12 );

 

save moham yexact1 yexact2 yapprox1 yapprox2 er1 er2 tt

B. Example 2

 

clear; clc; format long g; close all

 

T=1; dt=0.02; tt=0:dt:T; tic

 

%$\ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast$

 

h=0.2;

 

%$\ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast \ast$

 

u1exact=@(x)x;

 

u2exact=@(x)x;

 

ttt=0:h:T; alfa1=0; alfa2=0; v=6;

 

options = optimoptions('fsolve','Display', 'off',

 

'Algorithm', 'trust-region-dogleg',...

 

'MaxFunEvals',1e6,'MaxIter',1000,'TolFun',1e-10,'TolX',1e-10);

 

syms x t c1 c2 c3

 

u10=alfa1; u20=alfa2;

 

H1=c1+ c2$\ast$x+c3$\ast {\mathrm{x}}^{\wedge }$2;

 

h1=x-${\mathrm{x}}^{\wedge }$2; h2=x-(${\mathrm{x}}^{\wedge }$2/2)-($\mathrm{\hspace{0.17em}\hspace{0.17em}}{\mathrm{x}}^{\wedge }$3/3);

 

for i=1:length(ttt)-1

 

f10=subs(u10,x,t); f20=subs(u20,x,t); f30=f1$0^{\wedge }$2;

 

f1=f10+f20;

 

in1=int(f1,t,0,x );

 

u11=H1$\ast$( u10 -h1 -in1 );

 

f2=f30+f20;

 

in2=int(f2,t,0,x );

 

u21=H1$\ast$( u20 -h2 -in2 );

 

ut1=u10+u11; ut2=u20+u21;

 

g1=subs(ut1,x,t); g2=subs(ut2,x,t); g3=g$1^{\wedge }$2;

 

ig1=g1+g2; int1=int(ig1,t,0,x);

 

R1=ut1-h1-int1;

 

I1=vpa(R$1^{\wedge }$2,v);

 

I11=int(I1,x,ttt(i),ttt(i+1));

 

J1=vpa( I11, v);

 

fprintf('∖n∖n J1(x) = %s, ∖n',char((vpa(J1,v))) );

 

ig2=g2+g3; int2=int(ig2,t,0,x);

 

R2=ut2-h2-int2;

 

I2=vpa(R$2^{\wedge }$2,v);

 

I21=int(I2,x,ttt(i),ttt(i+1));

 

J2=vpa( I21, v);

 

fprintf('∖n∖n J2(x) = %s, ∖n',char((vpa(J2,v))) );

 

J11 = vpa(diff(J1,c1),v);

 

J12 = vpa(diff(J1,c2),v);

 

J13 = vpa(diff(J1,c3),v);

 

J21 = vpa(diff(J2,c1),v);

 

J22 = vpa(diff(J2,c2),v);

 

J23 = vpa(diff(J2,c3),v);

 

global dJ11 dJ12 dJ13 dJ21 dJ22 dJ23;

 

dJ11=@(c1,c2,c3)eval(char(J11)); dJ12=@(c1,c2,c3)eval(char(J12)); dJ13=@(c1,c2,c3)eval(char(J13));

 

dJ21=@(c1,c2,c3)eval(char(J21)); dJ22=@(c1,c2,c3)eval(char(J22)); dJ23=@(c1,c2,c3)eval(char(J23));

 

options = optimoptions('fsolve', 'Display', 'off',

 

'Algorithm', 'trust-region-dogleg',...

 

'MaxFunEvals',1e6,'MaxIter',1000,'TolFun',1e-10,'TolX',1e-10);

 

cc1=fsolve(@Sifun1,$[0\hspace{0.17em}\hspace{0.17em}0\hspace{0.17em}\hspace{0.17em}0]$,options);

 

cc2=fsolve(@Sifun2,$[0\hspace{0.17em}\hspace{0.17em}0\hspace{0.17em}\hspace{0.17em}0]$,options);

 

fprintf('∖n∖n∖n∖t∖t∖t∖t∖t∖tc11 = %1.10f,∖n∖n ∖t∖t∖t∖t∖t∖tc12 = %1.10f,∖n∖n ∖t∖t∖t∖t∖t∖tc13 = %1.10f,∖n∖n',cc1(1),cc1(2),cc1(3));

 

fprintf('∖n∖n∖n∖t∖t∖t∖t∖t∖tc21 = %1.10f,∖n∖n ∖t∖t∖t∖t∖t∖tc22 = %1.10f,∖n∖n ∖t∖t∖t∖t∖t∖tc23 = %1.10f, ∖n∖n',cc2(1),cc2(2),cc2(3));

 

u1=( subs( ut1,[c1,c2,c3],[cc1(1),cc1(2),cc1(3)] ) );

 

u2=( subs( ut2,[c1,c2,c3],[cc2(1),cc2(2),cc2(3)] ) );

 

fprintf('∖n∖n u1(%g) = %s, ∖n',ttt(i),char(vpa(u1,v)) );

 

fprintf('∖n∖n u2(%g) = %s, ∖n',ttt(i),char(vpa(u2,v)) );

 

yy1=@(x)eval(char(u1));

 

yy2=@(x)eval(char(u2));

 

yy=@(x)eval(char(u));

 

mt=ttt(i):dt:ttt(i+1);

 

for j=1:length(mt)

 

yyn1(j)=yy1(mt(j));

 

yyn2(j)=yy2(mt(j));

 

y1exn(j)=u1exact(mt(j));

 

y2exn(j)=u2exact(mt(j));

 

end

 

yapp$\{\mathrm{1},\mathrm{i}\}$=yyn1; yapp$\{\mathrm{2},\mathrm{i}\}$=yyn2; yex$\{\mathrm{1},\mathrm{i}\}$=y1exn; yex$\{\mathrm{2},\mathrm{i}\}$=y2exn;

 

u10=u1; u20=u2;

 

end

 

yapprox1=yapp$\{\mathrm{1,1}\}$; yapprox2=yapp$\{\mathrm{2,1}\}$; yexact1=yex$\{\mathrm{1,1}\}$; yexact2=yex$\{\mathrm{2,1}\}$;

 

for i=2:length(ttt)-1

 

yapprox1=[yapprox1,yapp$\{\mathrm{1},\mathrm{i}\}$(2:end)];

 

yapprox2=[yapprox2,yapp$\{\mathrm{2},\mathrm{i}\}$(2:end)];

 

yexact1=[yexact1,yex$\{\mathrm{1},\mathrm{i}\}$(2:end)];

 

yexact2=[yexact2,yex$\{\mathrm{2},\mathrm{i}\}$(2:end)];

 

end

 

th=0:0.1:T; er1=abs(yexact1-yapprox1); er2=abs(yexact2-yapprox2);

 

for i=1:length(th)

 

k(i)=round(th(i)/dt)+1;

 

end

 

L=$\left\{\hspace{0.17em}\right.$'t', 'Exact1', 'Exact2', ' MOHAM1', ' MOHAM2', 'Error1', 'Error2'$\left.\hspace{0.17em}\right\}$;

 

table(th', yexact1(k)', yexact2(k)', yapprox1(k)', yapprox2(k)',er1(k)',er2(k)', 'VariableNames',L)

 

figure('Units', 'characters', 'Name', 'p1', 'Position',

 

$[25\hspace{0.17em}\hspace{0.17em}6\hspace{0.17em}\hspace{0.17em}160\hspace{0.17em}\hspace{0.17em}45]$ );

 

plot(tt,yexact1,'linewidth',2);hold on;

 

plot(tt,yapprox1,'o','MarkerSize',8); grid on; grid minor; axis tight;

 

xlabel('x', 'fontweight', 'Bold');

 

legend('u_1 Exact solution', 'u_1 MOHAM solution', 'Location', 'northwest');

 

set( gca, 'fontname', 'Euclid', 'fontsize',12 );

 

figure('Units', 'characters', 'Name', 'p2', 'Position',

 

$[25\hspace{0.17em}\hspace{0.17em}6\hspace{0.17em}\hspace{0.17em}160\hspace{0.17em}\hspace{0.17em}45]$ );

 

plot(tt,yexact2,'linewidth',2);hold on;

 

plot(tt,yapprox2,'o','MarkerSize',8); grid on; grid minor; axis tight;

 

xlabel('x', 'fontweight', 'Bold');

 

legend('u_2 Exact solution', 'u_2 MOHAM solution', 'Location', 'northwest');

 

set( gca, 'fontname', 'Euclid', 'fontsize',12 );

 

figure('Units', 'characters', 'Name', 'Error1', 'Position',

 

$[25\hspace{0.17em}\hspace{0.17em}6\hspace{0.17em}\hspace{0.17em}160\hspace{0.17em}\hspace{0.17em}45]$ );

 

plot(tt,er1,'linewidth',2);

 

grid on; grid minor; axis tight; xlabel('x', 'fontweight', 'Bold');

 

ch=ylabel('u_1 Absolute Error', 'fontweight', 'Bold');

 

set( gca, 'fontname', 'Euclid', 'fontsize',12 );

 

figure('Units', 'characters', 'Name', 'Error2', 'Position',

 

$[25\hspace{0.17em}\hspace{0.17em}6\hspace{0.17em}\hspace{0.17em}160\hspace{0.17em}\hspace{0.17em}45]$ );

 

plot(tt,er2,'linewidth',2);

 

grid on; grid minor; axis tight; xlabel('x', 'fontweight', 'Bold');

 

ch=ylabel('u_1 Absolute Error', 'fontweight', 'Bold');

 

set( gca, 'fontname', 'Euclid', 'fontsize',12 );

 

save moham yexact1 yexact2 yapprox1 yapprox2 er1 er2 tt

References

[1] E. Gokmen, O. R. Isik, M. Sezer, "Taylor collocation approach for delayed Lotka-Volterra predator-prey system," Applied Mathematics and Computation, vol. 268 no. 1, pp. 671-684, DOI: 10.1016/j.amc.2015.06.110, 2015.

[2] K. Maleknejad, F. Mirzaee, "Numerical solution of linear Fredholm integral equations system by rationalized Haar functions method," International Journal of Computer Mathematics, vol. 80 no. 11, pp. 1397-1405, DOI: 10.1080/0020716031000148214, 2003.

[3] K. Maleknejad, F. Mirzaee, S. Abbasbandy, "Solving linear integro-differential equations system by using rationalized Haar functions method," Applied Mathematics and Computation, vol. 155 no. 2, pp. 317-328, DOI: 10.1016/S0096-3003(03)00778-1, 2004.

[4] Ş. Yüzbaşı, N. Şahin, M. Sezer, "A collocation approach to solving the model of pollution for a system of lakes," Mathematical and Computer Modelling, vol. 55 no. 3-4, pp. 330-341, DOI: 10.1016/j.mcm.2011.08.007, 2012.

[5] A. Arikoglu, I. Ozkol, "Solutions of integral and integro differential equation systems by using differential transform method," Computers & Mathematics with Applications, vol. 56 no. 9, pp. 2411-2417, DOI: 10.1016/j.camwa.2008.05.017, 2008.

[6] J. Biazar, E. Babolian, R. M. Islam, "Solution of a system of Volterra integral equations of the first kind by Adomian method," Applied Mathematics and Computation, vol. 139 no. 2-3, pp. 249-258, DOI: 10.1016/S0096-3003(02)00173-X, 2003.

[7] H. Aminikhah, J. Biazar, "A new analytical method for solving systems of Volterra integral equations," International Journal of Computer Mathematics, vol. 87 no. 5, pp. 1142-1157, DOI: 10.1080/00207160903128497, 2010.

[8] J. Biazar, M. Eslami, "Modified HPM for solving systems of Volterra integral equations of the second kind," Journal of King Saud University - Science, vol. 23 no. 1, pp. 35-39, DOI: 10.1016/j.jksus.2010.06.004, 2011.

[9] M. Javidi, "Modified homotopy perturbation method for solving system of linear Fredholm integral equations," Mathematical and Computer Modelling, vol. 50 no. 1-2, pp. 159-165, DOI: 10.1016/j.mcm.2009.02.003, 2009.

[10] J. Biazar, H. Ghazvini, "He’s homotopy perturbation method for solving systems of Volterra integral equations of the second kind," Chaos, Solitons & Fractals, vol. 39 no. 2, pp. 770-777, DOI: 10.1016/j.chaos.2007.01.108, 2009.

[11] R. Katani, S. Shahmorad, "Block by block method for the systems of nonlinear Volterra integral equations," Applied Mathematical Modelling, vol. 34 no. 2, pp. 400-406, DOI: 10.1016/j.apm.2009.04.013, 2010.

[12] K. Maleknejad, M. Tavassoli Kajani, "Solving linear integro-differential equation system by Galerkin methods with hydrid functions," Applied Mathematics and Computation, vol. 159 no. 3, pp. 603-612, DOI: 10.1016/j.amc.2003.10.046, 2004.

[13] K. Maleknejad, M. Shahrezaee, "Using Runge-Kutta method for numerical solution of the system of Volterra integral equation," Applied Mathematics and Computation, vol. 149 no. 2, pp. 399-410, DOI: 10.1016/S0096-3003(03)00148-6, 2004.

[14] M. Matinfar, M. Saeidy, J. Vahidi, "Application of homotopy analysis method for solving systems of Volterra integral equations," Advances in Applied Mathematics and Mechanics, vol. 4 no. 1, pp. 36-45, DOI: 10.4208/aamm.10-m1143, 2012.

[15] J. Pour-Mahmoud, M. Y. Rahimi-Ardabili, S. Shahmorad, "Numerical solution of the system of Fredholm integro-differential equations by the Tau method," Applied Mathematics and Computation, vol. 168 no. 1, pp. 465-478, DOI: 10.1016/j.amc.2004.09.026, 2005.

[16] J. Saberi-Nadjafi, M. Tamamgar, "The variational iteration method, a highly promising method for solving the system of integro differential equations," Computers & Mathematics with Applications, vol. 56 no. 2, pp. 346-351, DOI: 10.1016/j.camwa.2007.12.014, 2008.

[17] A. Akyuz-Dacoglu, M. Sezer, "Chebyshev polynomial solutions of systems of higher-order linear Fredholm–Volterra integro- differential equations," Journal of The Franklin Institute, vol. 342 no. 6, pp. 688-701, DOI: 10.1016/j.jfranklin.2005.04.001, 2005.

[18] E. Yusufoglu, "A homotopy perturbation algorithm to solve a system of Fredholm-Volterra type integral equations," Mathematical and Computer Modelling, vol. 47 no. 11-12, pp. 1099-1107, DOI: 10.1016/j.mcm.2007.06.022, 2008.

[19] N. R. Anakira, A. k. Alomari, A. F. Jameel, I. Hashim, "Multistage optimal homotopy asymptotic method for solving initial-value problems," Journal of Nonlinear Sciences and Applications, vol. 9 no. 4, pp. 1826-1843, DOI: 10.22436/jnsa.009.04.37, 2016.

[20] H. D. kasmaei, J. Rashidinia, "Optimal Homotopy Asymptotic and Homotopy Perturbation Methods for Linear Mixed Volterra-Fredholm Integral Equations," Nevşehir Bilim ve Teknoloji Dergisi, vol. 5 no. 2, pp. 86-103, DOI: 10.17100/nevbiltek.284733, 2016.

[21] M. S. Hashmi, N. Khan, S. Iqbal, "Numerical solutions of weakly singular Volterra integral equations using the optimal homotopy asymptotic method," Computers & Mathematics with Applications, vol. 64 no. 6, pp. 1567-1574, DOI: 10.1016/j.camwa.2011.12.084, 2012.

[22] V. Marinca, N. Herisanu, I. Nemes, "Optimal homotopy asymptotic method with application to thin film flow," Central European Journal of Physics, vol. 66 no. 3, pp. 648-653, DOI: 10.2478/s11534-008-0061-x, 2008.

[23] V. Marinca, N. Herisanu, C. Bota, B. Marinca, "An optimal homotopy asymptotic method applied to the steady flow of a fourth-grade fluid past a porous plate," Applied Mathematics Letters, vol. 22 no. 2, pp. 245-251, DOI: 10.1016/j.aml.2008.03.019, 2009.

[24] V. Marinca, N. Herişanu, "An optimal homotopy asymptotic method for solving nonlinear equations arising in heat transfer," International Communications in Heat and Mass Transfer, vol. 35 no. 6, pp. 710-715, DOI: 10.1016/j.icheatmasstransfer.2008.02.010, 2008.