\documentclass{mc}
%%%%% journal  info -- DO NOT CHANGE %%%%%%%%%
\setcounter{page}{197}
\renewcommand\thisnumber{x}
\renewcommand\thisyear {2021}
\renewcommand\thismonth{xxx}
\renewcommand\thisvolume{26}
\renewcommand\datereceived{August 27, 2020}
\renewcommand\dateaccepted{March 27, 2021}
%%%%%%%%   end journal info   %%%%%%%%%%%%%

\usepackage[displaymath,mathlines,pagewise]{lineno}

%Rad zaprimljen dana: 2020-08-27

%Rad prihvacen dana: 2021-03-27

%%%%% author macros %%%%%%%%%
% place your own macros HERE
%%%%% end %%%%%%%%%

%%%%%%%%%    already loaded packages -- DO NOT RELOAD  %%%%%
\usepackage[numbers,sort&compress]{natbib}
\usepackage{multirow}
\usepackage{multicol}
\usepackage{epstopdf}
\usepackage{enumerate}
\usepackage{float}
\usepackage[font=small,format=plain,labelfont=bf,up,textfont=normal,up]{caption}
\captionsetup[table]{skip=1pt}
\usepackage{soul}
\usepackage{caption}
\usepackage{subcaption}
\captionsetup[subtable]{font=footnotesize}
%%%%%%    end   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{document}
%\begin{linenumbers}
%%%%% title and author(s):
 \markboth{M.\,R.\,Doostdar, A.\,R.\,Vahidi, T.\,Damercheli and E.\,Babolian}{A hybrid functions method to Solve systems of ODEs}
% \title{Title}

%\markboth{D.\,Juki\'c}{How to use the {\sf mc.cls} class file}
\title{A hybrid functions method for solving linear and non-linear systems of ordinary differential equations}

% single author:
% \author[AUTHOR]{AUTHOR}
% \address{address of AUTHOR}
% \email{{\tt email address of AUTHOR} (AUTHOR)}

% multiple authors:
% Please mark \corrauth after the name of the corresponding author.
% different addresses:
%\author[AUTHOR1 and AUTHOR2]{AUTHOR1\affil{1}\comma\corrauth and AUTHOR2\affil{2}}
%\address{\affilnum{1}\ address of AUTHOR1\\   \affilnum{2}\ address of AUTHOR2}
%
%same address:
%Mohammadreza Doostdar, Alireza Vahidi, Tayebeh Damercheli, Esmail Babolian
\author[M.\,R.\,Doostdar, A.\,R.\,Vahidi, T.\,Damercheli, E.\,Babolian]{Mohammadreza Doostdar\affil{1}, Alireza Vahidi\affil{1}\corrauth,\linebreak Tayebeh Damercheli\affil{1} and Esmail Babolian\affil{2}}
\address{\affilnum{1}{Department of Mathematics, College of Science, Yadegar-e-Imam Khomeyni (RAH) Shahr-e-Rey Branch, Islamic Azad University, Persian Gulf Freeway, Tehran, 1815163111, Iran}\\ \affilnum{2}{Faculty of Mathematical Sciences and Computer, Kharazmi University, 50 Taleghani Avenue, Tehran, 1561836314, Iran}}
%
\emails{{\tt mrdoostdar@yahoo.com}\,\,(M.\,R.\,Doostdar), {\tt alrevahidi@yahoo.com} (A.\,R.\,Vahidi), {\tt tdamercheli@gmail.com} (T.\,Damercheli),\newline {\tt babolian@khu.ac.ir} (E.\,Babolian)}
%\emails{{\tt email of AUTHOR1} (AUTHOR1), {\tt email of AUTHOR2} (AUTHOR2), {\tt email of AUTHOR3} (AUTHOR3)}

%%%%% Begin Abstract %%%%%%%%%%%
\begin{abstract}
In the present paper, we use a hybrid method to solve linear or non-linear systems of ordinary differential equations (ODEs). By using this method, these systems are reduced to a linear or non-linear system of algebraic equations. In error discussion of the suggested method, an upper bound of the error is obtained. Also, to survey the accuracy and the efficiency of the present method, some examples are solved and comparisons between the obtained results with those of several other methods are carried out.
\end{abstract}
%%%%% end %%%%%%%%%%%

%%%%% Keywords %%%%%%%%%%%
\keywords{block-pulse functions, hybrid functions, Legendre polynomials, system of ordinary differential equations}

%%%% AMS subject classifications %%%%
\ams{34A30, 34A34}

%%%% maketitle %%%%%
\maketitle

%%%% Start %%%%%%
\section{Introduction}\label{sec:1}
Studies on differential equations and their systems have always been considered due to their widespread applications in physical systems, engineering sciences, medicine, economy, biology, population systems,  etc. \cite{AbellBraselton,Butcher}. For instance, the predator and prey problem in biological systems  \cite{BiazarMontazeri2005, Dubeyetal2001,Lotka}, the mathematical models for describing   CD4$^+$ T cells and HIV interactions in medicine \cite{culshaw2000,essunger1994, kirschner1996,nowak1996,nowakl1991,perelson1990,perelson1993,lwang2006}, and stiff problems in biochemistry and life sciences \cite{FieldNoyes1974}, are expressed as linear or non-linear systems of ODEs. Many researchers have applied various methods for solving systems of ODEs. The Adomian decomposition method (ADM) \cite{Biazaretal2004,Mahmoodetal2005}, the variational iteration method (VIM) \cite{Biazar2007,Darvishietal2007,merdan2011}, the homotopy perturbation method (HPM) \cite{Biazar2015,merdan2007}, the differential transform method \cite{srivastava2014,ThongmoonPusjuso2010}, and a collocation method \cite{Rahimietal2017, yuzbasi2012, yuzbasi2016} have been proposed for solving linear and non-linear systems of ODEs.

In recent years, using hybrid functions has been considered for solving various mathematical models \cite{Asadyetal2005,HesameddiniShahbazi2018,Marzbanetal2009,MarzbanRazzaghi2006,Marzbanetal2008,MashayekhiRazzaghi2016,Mirzaeeetal2019,Mohammadietal2017,TavassoliVencheh2005}. One of these functions is obtained by the combination of block-pulse functions (BPFs) and Legendre polynomials. Based on a hybrid of BPFs and Legendre polynomials, Hsiao \cite{hsiao2009} solved the  Fredholm and Volterra integral equations; Singh et al. \cite{Singh2010} introduced a numerical method for numerical evaluation of the Hankel transform; a Maleknejad and Hashemizadeh \cite{MaleknejadHashemi2011} presented an approach to solve Hammerstein integral equations of mixed type; Maleknejad and Ebrahimzadeh \cite{MaleknejadEbrahimzadeh2014} proposed a method for solving optimal control of Volterra integral systems. Also, by using these hybrid functions Sahu and Saha Ray \cite{Sahu2015} solved a system of non-linear Fredholm-Hammerstein; Rafiei et al. \cite{Rafieietal2017} found the optimal solution of linear time delay systems; Maleknejad and Saeedipoor \cite{MaleknejadSaeedipoor2017} solved   Fredholm integral equation of the first kind; Hesameddini and Riahi \cite{Hesameddini2019} obtained a numerical solution of partial differential equations with non-local integral conditions, and finally Saha Ray and Singh \cite{SahaSingh2018} found  the numerical solutions of stochastic Volterra-Fredholm integral equations.
Here, considering the general form of the system of ODEs as
\begin{equation}\label{system}
\left \{
\begin{aligned}
  \frac{\mathrm{d}u_1}{\mathrm{d}t}&=f_1(t, u_1, u_2, \ldots, u_n),\\
   \frac{\mathrm{d}u_2}{\mathrm{d}t}&=f_2(t, u_1, u_2, \ldots, u_n),\\
 &\ \vdots\\
 \frac{\mathrm{d}u_n}{\mathrm{d}t}&=f_n(t, u_1, u_2, \ldots, u_n),
\end{aligned}
\right.
\end{equation}
with the initial conditions $u_1(0)=\lambda_1$, $u_2(0)=\lambda_2$, \ldots, $u_n(0)=\lambda_n$, we apply the hybrid BPFs and Legendre polynomials method (HBPLM) to solve this system. By using this method, systems of ODEs are reduced to a linear or non-linear system of algebraic equations, which can be solved by the proper methods. In addition, since every ordinary differential equation can be converted into a system of differential equations of the first order, the HBPLM can be implemented to solve most differential equations of higher orders.

We arranged the remainder of the paper as follows: In Section \ref{sec:2}, we review required definitions and some properties of hybrid BPFs and Legendre polynomials.  In Section \ref{sec:3}, we apply the HBPLM for linear and non-linear systems of ODEs. Error analysis of the presented method is studied  in Section \ref{sec:4}. Finally, in Section \ref{sec:5}, to verify the effectiveness of the proposed method, a mathematical model for HIV interactions with CD4$^+$ T cells, the predator and prey problem, some examples of stiff systems, and two other examples will be solved and the obtained absolute errors are compared with those of several other methods.
\section{Preliminaries}\label{sec:2}
In this section, required definitions and some properties of hybrid BPFs and Legendre polynomials and operational matrices are reviewed.
\subsection{Hybrid BPFs and Legendre polynomials}%subsection2-1 ------------------------------------------------------------------------------------------------
Considering the interval $[0, T_f)$ and by   hybridization of BPFs and the well-known Legendre polynomials denoted by $L_q(t)$, on this interval, hybrid functions $\ell_{ij}, i=1, 2, \ldots, N$ and $j=0, 1, \ldots, M-1$ are defined as \cite{hsiao2009}
%\begin{equation}
\[
\ell_{ij}(t)=\left\{
\begin{aligned}
 &L_j(\frac{2N}{T_f}t-2i+1), & t\in I_i,\\
 &0, & \text{o. w,}
\end{aligned}
\right.
\]
%\end{equation}
where $N$ and $M$ indicate the order of BPFs and the order of Legendre polynomials, respectively,  and $I_i=[\frac{i-1}{N}T_f,\frac{i}{N}T_f)$. Without loss of generality, we let $T_f=1$.\\
Any function $f\in L^2[0, 1)$ can be expanded by the basis functions  $\{\ell_{ij}\}$ as
\begin{equation}\label{1}
f(t)=\sum_{i=1}^{\infty}\sum_{j=0}^{\infty}\alpha_{ij}\ell_{ij}(t),
\end{equation}
where $\alpha_{ij}$ are   hybrid coefficients.
By truncating the infinite series in equation \eqref{1} at some values of $N$ and $M$, it can be written as
\begin{equation}\label{fnm}
f(t)\simeq f_{NM}(t)=\sum_{i=1}^{N}\sum_{j=0}^{M-1}\alpha_{ij}\ell_{ij}(t)=\textbf{A}^T\textbf{L}(t),
\end{equation}
where \textbf{A} and \textbf{L}(t) are defined as follows:
%\begin{equation}
\[
\textbf{A}=[\alpha_{10}, \alpha_{11}, \dots, \alpha_{1(M-1)}, \alpha_{20}, \alpha_{21}, \dots, \alpha_{2(M-1)}, \dots, \alpha_{N0}, \alpha_{N1}, \dots, \alpha_{N(M-1)}]^T,\label{A}
\]
%\end{equation}
and
%\begin{equation}
\[
\textbf{L}(t)=[\ell_{10}(t),\,\dots, \ell_{1(M-1)}(t), \ell_{20}(t),\,\dots, \ell_{2(M-1)}(t), \dots, \ell_{N0}(t),\,\dots, \ell_{N(M-1)}(t)]^T.
\]
%\end{equation}
\subsection{Operational matrices} % subsection 2-2-----------------------------------------------------------------------------------------------------------------------------
Suppose that $P$ denotes the operational matrix of the integration of the vector \textbf{L}(t) defined above; then
%\begin{align}
\[
\int_0^t \textbf{L}(\eta)\, \mathrm{d} \eta \simeq P\textbf{L}(t),
\]%\end{align}
where $P$ is an $NM\times NM$ matrix   defined as \cite{hsiao2009}:
%\begin{equation}
\[P =
\begin{pmatrix}
D & S & S & \ldots & S\\
\textbf{0} & D & S & \ldots & S \\
\textbf{0} &\textbf{0} & D & \ldots & S \\
\vdots & \vdots & \vdots & \ddots & \vdots\\
\textbf{0}&\textbf{0}&\textbf{0}&\ldots&D
\end{pmatrix}
,
\]%\end{equation}
where the zero matrix $\textbf{0}$, $S$ and $D$ are $M\times M$ matrices and
%\begin{equation}
\[S=\frac{T_f}{N}
\begin{pmatrix}
1 & 0 & 0 & \ldots & 0\\
0 & 0 & 0 & \ldots & 0 \\
0 & 0 & 0 & \ldots & 0 \\
\vdots & \vdots & \vdots & \ddots & \vdots\\
0 & 0 & 0 &\ldots & 0
\end{pmatrix}
,
\]%\end{equation}
\[%\begin{equation}
D =\frac{T_f}{2N} \left (\begin{array}{*{11}{c}}
1 & 1 & 0 & 0 & 0 & \ldots & 0 & 0 & 0 & 0 & 0 \\
-\frac{1}{3} & 0 & \frac{1}{3} & 0 & 0 & \ldots & 0 & 0 & 0 & 0 & 0\\
0 & -\frac{1}{5} & 0 & \frac{1}{5} & 0 & \ldots & 0 & 0 & 0 & 0 & 0\\
0 & 0 & -\frac{1}{7} & 0 & \frac{1}{7}& \ldots & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & -\frac{1}{9} & 0 & \ldots & 0 & 0 & 0 & 0 & 0 \\
\vdots & \vdots & \vdots & \vdots &  \vdots& \ddots & \vdots & \vdots & \vdots & \vdots & \vdots \\
0 & 0 & 0 & 0 & 0 & \ldots & 0 & \frac{1}{2M-9} & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & \ldots & \frac{-1}{2M-7} & 0 & \frac{1}{2M-7} & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & \ldots & 0 &  \frac{-1}{2M-5} & 0 &  \frac{1}{2M-5} & 0 \\
0 & 0 & 0 & 0 & 0 & \ldots & 0 & 0 &  \frac{-1}{2M-3} & 0&  \frac{1}{2M-3}\\
0 & 0 & 0 & 0 & 0 & \ldots & 0 & 0 & 0 & \frac{-1}{2M-1} & 0
\end{array}\right ).
\]%\end{equation}
To evaluate the product of $\textbf{L}(t)$ and $\textbf{L}^T(t)$, let
\[%\begin{equation}\label{3}
\textbf{L}(t) \textbf{L}^T(t)\textbf{A}=\tilde{\textbf{A}}\textbf{L}(t),
\]%\end{equation}
where $\tilde{\textbf{A}}$ is an $NM\times NM$ matrix   defined in \cite{hsiao2009}.
%----------------------------------------------------------------------------------------------------------------------------------------------------------------------
\section{The HBPLM for linear and non-linear systems of ODEs}\label{sec:3}
In this section, the HBPLM will be implemented to approximate the solutions of linear and non-linear systems of ODEs.
\subsection{Implementation for linear systems of ODEs}%-------------------------------------------Subsection3-1
We consider linear form of   system \eqref{system} as
\begin{equation}\label{systemlinear}
\left \{
\begin{aligned}
 \frac{\mathrm{d}u_1}{\mathrm{d}t}&=g_1(t)+\sum_{k=1}^nh_{1,k}(t)u_k(t),\\
 \frac{\mathrm{d}u_2}{\mathrm{d}t}&=g_2(t)+\sum_{k=1}^nh_{2,k}(t)u_k(t),\\
&\ \vdots\\
 \frac{\mathrm{d}u_n}{\mathrm{d}t}&=g_n(t)+\sum_{k=1}^nh_{n,k}(t)u_k(t),
\end{aligned}
\right.
\end{equation}
where $g_k$, $h_{r,k}$  for $k,r=1, 2, \ldots, n$, are the known functions.
From equation \eqref{fnm}, $u_k(t)$ can be written as
\begin{equation}\label{expansionx}
u_k(t)\simeq \textbf{A}_k^T\textbf{L}(t),\qquad k=1, 2, \ldots, n,
\end{equation}
where the unknown vector $\textbf{A}_k$, which must be calculated, is considered as
\begin{equation}
\begin{aligned}
&\textbf{A}_k=[a_{k,1}\quad a_{k,2}\quad \ldots \quad a_{k,NM}]^T.
\end{aligned}
\end{equation}
Now, we approximate $u_k'(t)$ as follows:
\begin{equation}\label{expansionxprim}
u_k'(t)\simeq {\textbf{A}}_{k}'^T\textbf{L}(t), \qquad k=1, 2, \ldots, n,
\end{equation}
where ${\textbf{A}}_{k}'\simeq(\textbf{A}_{k}^T-u_{k(0)}^T)P^{-1}$ and $u_{k(0)}$ can be obtained as follows:
\begin{equation}
u_{k(0)}=\bigg(\begin{matrix} \underbrace{\begin{matrix}u_k(0)&\overbrace{\begin{matrix}0& 0&\ldots& 0\end{matrix}}^{M-1}& u_k(0)&\overbrace{\begin{matrix}0& 0&\ldots& 0\end{matrix}}^{M-1}& \ldots& u_k(0)&\overbrace{\begin{matrix}0& 0&\ldots& 0\end{matrix}}^{M-1}\end{matrix}}_{NM}\end{matrix}\bigg)^T.
\end{equation}
Considering equation \eqref{fnm}, the known functions $g_k$ and $h_{r,k},\, r,k=1, 2, \ldots, n$, by the basis functions $\{\ell_{ij}\}$, can be written as
\begin{equation}\label{expansiong}
g_k(t)\simeq \mathcal{G}_k^T\textbf{L}(t),
\end{equation}
and
\begin{equation}\label{app_h}
h_{r,k}(t)\simeq \mathcal{H}_{r,k}^T\textbf{L}(t).
\end{equation}
Also, the terms $h_{r,k}(t)u_k(t),\, r,k=1, 2, \ldots, n$, can be expanded in terms of the vector $\textbf{L}(t)$ as
\begin{equation}\label{htx}
\begin{aligned}
h_{r,k}(t)u_k(t)\simeq (\mathcal{H}_{r,k}^T\textbf{L}(t))(\textbf{A}_{k}^T\textbf{L}(t))=\mathcal{H}_{r,k}^T\textbf{L}(t)\textbf{L}^T(t)\textbf{A}_{k}
=\mathcal{H}_{r,k}^T\tilde{\textbf{A}}_{k}\textbf{L}(t).
\end{aligned}
\end{equation}
Note that if $h_{r,k}(t)$ is  a constant as $c$, then equation \eqref{htx} can be written as
\begin{equation*}
h_{r,k}(t)u_k(t)\simeq c\textbf{A}_k^T\textbf{L}(t), \qquad r, k=1, 2, \ldots, n.
\end{equation*}
Now, by substituting equations \eqref{expansionx} - \eqref{htx} into equation \eqref{systemlinear} and replacing $\simeq$ with $=$, we have
\begin{equation}\label{systemlinear-imp}
\left \{
\begin{aligned}
 \bigg((\textbf{A}_{1}^T-u_{1(0)}^T)P^{-1}-\mathcal{G}_1^T-\sum_{k=1}^n\mathcal{H}_{1,k}^T\tilde{\textbf{A}}_{k}\bigg)\textbf{L}(t)&=0,\\
 \bigg((\textbf{A}_{2}^T-u_{2(0)}^T)P^{-1}-\mathcal{G}_2^T-\sum_{k=1}^n\mathcal{H}_{2,k}^T\tilde{\textbf{A}}_{k}\bigg)\textbf{L}(t)&=0,\\
&\ \vdots\\
 \bigg((\textbf{A}_{n}^T-u_{n(0)}^T)P^{-1}-\mathcal{G}_n^T-\sum_{k=1}^n\mathcal{H}_{n,k}^T\tilde{\textbf{A}}_{k}\bigg)\textbf{L}(t)&=0.
\end{aligned}
\right.
\end{equation}
By collocating system \eqref{systemlinear-imp} at the points $t=\frac{2\zeta-1}{2NM}$, $\zeta=1,2, \ldots, NM$, and solving the obtained linear system of algebraic equations by a well-known numerical method, the vector $\textbf{A}_k$ is obtained and $u_k(t)$ is consequently determined by using equation \eqref{expansionx}.
\subsection{Implementation for non-linear systems of ODEs}%-------------------------------------------Subsection3-2
A non-linear system of ODEs is considered as
\begin{equation}\label{systemnonlinear}
\left \{
\begin{aligned}
 \frac{\mathrm{d}u_1}{\mathrm{d}t}&=f_{1}(t, u_{1}, u_{2}, \ldots, u_n),\\
 \frac{\mathrm{d}u_2}{\mathrm{d}t}&=f_2(t, u_1, u_2, \ldots, u_n),\\
&\ \vdots\\
 \frac{\mathrm{d}u_n}{\mathrm{d}t}&=f_n(t, u_1, u_2, \ldots, u_n).
\end{aligned}
\right.
\end{equation}
Here, $u_k(t)$ and $u'_k(t)$ are expanded as equations  \eqref{expansionx} and \eqref{expansionxprim}, respectively. If necessary, we use equation \eqref{expansiong} to expand a known function. Also, non-linear terms can be expanded in terms of the vector $\textbf{L}(t)$ as
\[%\begin{equation}\label{TI}
\begin{aligned}
u_r(t)u_k(t) \simeq (\textbf{A}_{r}^T\textbf{L}(t))(\textbf{A}_{k}^T\textbf{L}(t))=\textbf{A}_{r}^T\textbf{L}(t)\textbf{L}^T(t)\textbf{A}_{k}=\textbf{A}_{r}^T\tilde{\textbf{A}}_{k}\textbf{L}(t),\, r,k=1, 2, \ldots, n,
\end{aligned}
\]%\end{equation}
and
\[%\begin{equation}\label{T3}
\begin{aligned}
u_k^3(t) =u_k^2(t).u_k(t)\simeq (\textbf{A}_{k}^T\tilde{\textbf{A}}_{k}\textbf{L}(t))(\textbf{A}_{k}^T\textbf{L}(t))=\textbf{A}_{k}^T\tilde{\textbf{A}}_{k}\textbf{L}(t)\textbf{L}^T(t)\textbf{A}_{k}
=\textbf{A}_{k}^T\tilde{\textbf{A}}_{k}^2\textbf{L}(t),
\end{aligned}
\]%\end{equation}
and in the same way, $u_k^{\alpha}(t) \simeq \textbf{A}_{k}^T\tilde{\textbf{A}}_{k}^{\alpha-1}\textbf{L}(t)$, where $\alpha$ is a positive integer. Suppose
\[%\begin{equation}\label{KLt}
f_k(t, u_1, u_2, \ldots, u_n)\simeq \bar{\textbf{F}}_k^T\textbf{L}(t), \qquad k=1, 2, \ldots, n;
\]%\end{equation}
then  system \eqref{systemnonlinear} can be written as the following non-linear algebraic system:
\begin{equation}\label{system3}
\bigg((\textbf{A}_k^T-u_{k(0)}^T)P^{-1}-\bar{\textbf{F}}_k^T\bigg)\textbf{L}(t)=0, \qquad k=1, 2, \ldots, n.
\end{equation}
Like the linear form, by collocating system \eqref{system3} at the points $t=\frac{2\zeta-1}{2NM}$, $\zeta=1,2, \ldots, NM$, and solving the obtained  non-linear system of algebraic equations by a proper method such as Newton's method, the vector $\textbf{A}_k$ is obtained and $u_k(t)$ is consequently determined by using equation \eqref{expansionx}.
%------------------------------------------------End of Section{Implementation of the method}
\section{Error analysis}\label{sec:4}
In this section, considering the Sobolev space and the associated norm, an upper bound of the error for the HBPLM is obtained. We consider the Sobolev norm on the interval $(-1,1)$ as follows:
\begin{equation*}
\|f\|_{H^d(-1,1)}=\bigg(\sum_{i=0}^d\|f^{(i)}\|_{L^2(-1,1)}^2\bigg)^{\frac{1}{2}},
\end{equation*}
where $d\ge 0$ is an integer.
\begin{lemma}\label{lem1}\textup{\cite{canuto.etal}}
Let $\{L_q\}_{q=0}$ be the sequence of Legendre polynomials and $P_S(t)=\sum_{q=0}^Sa_qL_q(t)$   the best polynomial approximation of degree $S$ for $f\in L^2(-1,1)$. Then, for $d\ge 1$, there exists a constant $C_0>0$ such that
\[%\begin{equation}
\|f-P_S\|_{L^{\infty}(-1,1)}\leq C_0S^{\frac{3}{4}-d}\|f\|_{H^d(-1,1)},
\]%\end{equation}
for all functions $f$ in $H^d(-1,1)$.
\end{lemma}
\begin{lemma}\label{lem2}
Let $f\in H^d[0,1)$ and let $f_{NM}$ be the polynomial approximation of $f$ defined in equation \eqref{fnm}. Then,
\[%\begin{equation}
\|f-f_{NM}\|_{L^{\infty}[0,1)}\leq C_0(NM)^{\frac{3}{4}-d}\max\limits_{1\leq i\leq N}\|f\|_{H^d(I_i)},
\]%\end{equation}
where $I_i=[\frac{i-1}{N},\frac{i}{N})$ and $C_0$ is a positive constant.
\end{lemma}
\begin{proof}
It can be obviously concluded by using Lemma \ref{lem1}.
\end{proof}
\noindent Before we express the following theorem, we represent system \eqref{system} as
\begin{equation}\label{intequation}
U(t)=U(0)+\int_0^t F(\tau, U(\tau))\, \mathrm{d} \tau,
\end{equation}
where
\begin{align*}
 U(t)&=[u_1(t), u_2(t),\ldots, u_n(t)]^T,\\
 F(\tau, U(\tau))&=
\begin{pmatrix}
f_1(\tau, U(\tau))\\
f_2(\tau, U(\tau))\\
\vdots\\
f_n(\tau, U(\tau))
\end{pmatrix}.
\end{align*}
\begin{theorem}
Let $U\in H^d[0,1)$ be the exact solution of equation \eqref{intequation} and $\tilde{U}(t)=U_{NM}(t)$   the approximate solution obtained by the HBPLM. Moreover, assume that $f_k(\tau, U(\tau)),\, k=1, 2, \ldots, n$ is a continuous function for $0\leq \tau\leq t<1$ and that satisfies the Lipschitz condition
\begin{equation}\label{Lipcon}
|f_k(\tau, U(\tau))-f_k(\tau, W(\tau))|\leq L_k\,\|U-W\|_{\infty},
\end{equation}
where $L_k>0$, $k=1, 2, \ldots, n$, is a Lipschitz constant. Then, there exists a constant $\delta>0$ such that
\[%\begin{equation}
\|U-\tilde{U}\|_{\infty}\leq\delta \,\max\limits_{1\leq i\leq N}\|u_p\|_{H^d(I_i)}.
\]%\end{equation}
\end{theorem}
\begin{proof}
Let $\tilde{u}_k(t)=u_{k(NM)}(t)$, $k=1, 2, \ldots, n$, be the approximate solution of the system defined in equation \eqref{intequation}, and let  $e_k(t)=u_k(t)-\tilde{u}_k(t)$ be the error term. Then,
\[%\begin{equation}\label{error}
e_k(t)=u_k(t)-\tilde{u}_k(t)=\int_0^t \big(f_k(\tau, U(\tau))-f_k(\tau, \tilde{U}(\tau))\big)\,\mathrm{d}\tau.
\]%\end{equation}
Using equation \eqref{Lipcon} and $0\leq t<1$, we obtain
\[%\begin{equation}\label{errornorm}
\begin{aligned}
|e_k(t)|\leq L_k\,\|U-\tilde{U}\|_{\infty}=L_k\max_{\substack{1\leq k\leq n}}\|u_k-\tilde{u}_k\|_{\infty}.
\end{aligned}
\]%\end{equation}
Let $L=\max\limits_{1\leq k\leq n}\,L_k$; then
\begin{equation}\label{42}
\|e_k\|_{\infty}\leq L\,\|u_p-\tilde{u}_p\|_{\infty},\qquad p\in\{1, 2, \ldots, n\}.
\end{equation}
Now, by using Lemma \ref{lem2} and equation \eqref{42}, we have
\begin{equation}\label{43}
\|e_k\|_{\infty}\leq L\,\|u_p-\tilde{u}_p\|_{\infty}=L\,\|u_p-\tilde{u}_p\|_{L^{\infty}[0,1)}\leq \delta\,\max\limits_{1\leq i\leq N}\|u_p\|_{H^d(I_i)},
\end{equation}
where $\delta=L C_0 (NM)^{\frac{3}{4}-d}$.\\
If $e(t)=U(t)-\tilde{U}(t)$, then using \eqref{43}, we get
\[%\begin{equation}
\|e\|_{\infty}=\|U-\tilde{U}\|_{\infty}\leq \delta\,\max\limits_{1\leq i\leq N}\|u_p\|_{H^d(I_i)}.
\]%\end{equation}
\end{proof}
%------------------------------------------------% Section Numerical Example
\section{Numerical examples}\label{sec:5}
\noindent In this section, the effectiveness of the present method is studied by applying the HBPLM to several examples of linear and non-linear systems of ODEs. As the first example, a mathematical model for HIV interactions with CD4$^+$ T cells is considered. Then, the predator-prey problem, two examples of linear and non-linear stiff systems, and two other examples are solved. All computations are performed using Matlab 2017a software package on a laptop with the Intel core i5-3210M CPU processor and 4GB RAM.
%----------------Example 1
\begin{example}
\emph{As the first example, we consider a mathematical model of HIV interactions with CD4$^+$ T cells, which is described by a three-dimensional system of non-linear ODEs as \cite{merdan2011,ongun2011,srivastava2014,lwang2006,yuzbasi2012}:
\begin{equation}\label{systemHIV}
\left \{
\begin{aligned}
T'(t)&=s-\mu T+rT\bigg(1-\frac{T+I}{T_{max}}\bigg)-\alpha VT,\\
 I'(t)&=\alpha VT-\beta I,\\
 V'(t)&=C\beta I- \gamma V,
\end{aligned}
\right.
\end{equation}
where the initial conditions and parameters are as given
\begin{equation*}
\left \{
\begin{aligned}
 T(0)&=0.1,  \quad &I(0)&=0, \quad &V(0)&=0.1,\\
 r&=3, \quad &s&=0.1, \quad &\mu&=0.02, \quad &T_{max}&=1500, \quad &\alpha&=0.0027,\\
 \beta&=0.3, \quad &C&=10, \quad &\gamma&=2.4.
\end{aligned}
\right.
\end{equation*}
In this example, the HBPLM is implemented for the above model for $N=2,M=8$, which reduces it to a non-linear system of algebraic equations. Since system \eqref{systemHIV} has no exact solution, we prefer to use the solution of the classical forth-order Runge-Kutta method (RK4) as an acceptable solution for the purpose of comparison. Hence,   numerical results of the suggested method and those of the modified variational iteration method (MVIM) \cite{merdan2011}, the Bessel collocation method (BCM) \cite{yuzbasi2012}, the differential transform method (DTM) \cite{srivastava2014}, and the Laplace Adomian decomposition method (LADM) \cite{ongun2011} for $T(t)$, $I(t)$ and $V(t)$, are compared with the RK4  solutions. Tables \ref{tb1:subtable1}, \ref{tb2:subtable2} and \ref{tb3:subtable3} show these comparisons.
% Tables example 1---------------------
\begin{table}[h]
\centering
\begin{subtable}{.75\textwidth}
\centering
\subcaption{Comparisons for $T(t)$}
\begin{tabular}{cccccc}
\hline
\multirow{1}{*}{t} & \multirow{1}{*}{HBPLM} &  \multirow{1}{*}{MVIM \cite{merdan2011}}&   \multirow{1}{*}{LADM \cite{ongun2011}}&
 \multirow{1}{*}{BCM \cite{yuzbasi2012}}& \multirow{1}{*}{DTM \cite{srivastava2014}}\\
\hline
0.0 & 0 & 0 & 0 & 0 & 0 \\
0.2 & 1.80E-09 & 3.50E-09&8.10E-07 & 4.95E-03& 2.84E-03 \\
0.4 & 0.90E-09 & 2.56E-07&1.35E-04&2.59E-02&1.64E-02 \\
0.6 & 0.90E-09&0.00E-00&3.28E-03&6.90E-02&5.35E-02\\
0.8 & 21.9E-09 & 0.00E-00& 3.67E-02& 1.38E-01& 1.32E-01\\
\hline
\end{tabular}
\label{tb1:subtable1}
\end{subtable}\
\begin{subtable}{.75\textwidth}
\centering
\subcaption{Comparisons for $I(t)$}\
\begin{tabular}{cccccc}
\hline
\multirow{1}{*}{t} & \multirow{1}{*}{HBPLM} &  \multirow{1}{*}{MVIM \cite{merdan2011}}&   \multirow{1}{*}{LADM \cite{ongun2011}}&
 \multirow{1}{*}{BCM \cite{yuzbasi2012}}& \multirow{1}{*}{DTM \cite{srivastava2014}}\\
\hline
0.0 & 0 & 0 & 0 & 0 & 0 \\
0.2 & 9.10E-14 & 4.99E-13 & 5.13E-12 & 2.15E-07 & 3.34E-07 \\
0.4 & 1.90E-13 & 3.91E-11 & 8.21E-10 & 2.23E-07 & 8.34E-07 \\
0.6 & 3.70E-13 & 4.75E-10 & 4.46E-08 & 8.71E-07 & 1.43E-06\\
0.8 & 5.60E-13 & 2.91E-09 & 1.08E-07 & 1.80E-06 & 3.11E-06\\
\hline
\end{tabular}
\label{tb2:subtable2}
\end{subtable}
\begin{subtable}{.75\textwidth}
\centering
\subcaption{Comparisons for $V(t)$}
\begin{tabular}{cccccc}
\hline
\multirow{1}{*}{t} & \multirow{1}{*}{HBPLM} &  \multirow{1}{*}{MVIM \cite{merdan2011}}&   \multirow{1}{*}{LADM \cite{ongun2011}}&
 \multirow{1}{*}{BCM \cite{yuzbasi2012}}& \multirow{1}{*}{DTM \cite{srivastava2014}}\\
\hline
0.0 & 0 & 0 & 0 & 0 & 0 \\
0.2 & 1.00E-10 & 6.54E-08 & 1.17E-07 & 7.52E-08 & 1.57E-07 \\
0.4 & 2.00E-10 & 1.07E-06 & 1.84E-05 & 4.71E-08 & 1.41E-05 \\
0.6 & 1.00E-10 & 5.74E-06 & 6.87E-04 & 2.31E-07 & 2.15E-04\\
0.8 & 1.00E-10 & 2.01E-05 & 4.71E-03 & 7.94E-07 & 1.53E-03\\
\hline
\end{tabular}
\label{tb3:subtable3}
\end{subtable}
\vspace*{0.2cm}
\caption{\it Comparisons between  approximate solutions for (a) $T(t)$, (b) $I(t)$ and (c) $V(t)$ by \emph{HBPLM, MVIM, LADM, BCM} and \emph{DTM} with \emph{RK4}}\label{TableT}
\end{table}
}
\end{example}

%*********** Example 2
\begin{example}\label{EX2}
\emph{Consider the following system of ODEs that arises in the modeling of the two species predator and prey problem \cite{Butcher}:
\begin{equation}\label{example2}
\left \{
\begin{aligned}
 u'(t)&=\big(2-v(t)\big)u(t),\\
 v'(t)&=\big(u(t)-1\big)v(t),\\
\end{aligned}
\right.
\end{equation}
in which $u$ denotes the prey population and $v$ is the population of  predators. With the passage of time, a permanently repeated cycle of interrelated falls and rises in the predator and prey populations occurs. The time period of these repeated cycles, which is indicated by the letter $\mathcal{T}$, has been calculated in \cite{Butcher}, from which, with starting values $u(0)=2$ and $v(0)=2$, we see that the time period converges to $\mathcal{T}\simeq4.61$. We solved system \eqref{example2} by the HBPLM for $N=2, M=8$ and the  RK4 method. Figure \ref{Fig1} shows the approximate solutions obtained by the HBPLM and RK4. Also, the obtained numerical results are shown in Table \ref{Table-ex2}.
\begin{table}[H]
\centering
\resizebox{1\textwidth}{!}{
\begin{tabular}{lccc|ccc}
\hline
\multirow{2}{*}{t} & \multicolumn{3}{c}{$u$} &  \multicolumn{3}{c}{$v$}\\
\cline{2-4}
\cline{5-7}
& RK4 & HBPLM & absolute error&RK4 & HBPLM & absolute error \\
\hline
0.5 & 1.531714 & 1.529849 & 1.87E-03 & 3.036906 &3.038620  & 1.71E-03 \\
1.0 & 0.819587 & 0.819081 & 5.06E-04 & 3.274107 & 3.275590 & 1.48E-03 \\
1.5 & 0.490400 & 0.491720 & 1.32E-03 & 2.715792 & 2.713292 & 2.50E-03\\
2.0 & 0.406754 & 0.404548 & 2.21E-03 & 2.047384 & 2.050809 & 3.43E-03\\
2.5 & 0.455341 & 0.449063 & 6.28E-03 & 1.533120 & 1.542198 & 9.08E-03\\
3.0 & 0.628130 & 0.627862 & 2.67E-04 & 1.212053 & 1.211689 & 3.64E-04\\
3.5 & 0.969317 & 0.969139 & 1.78E-04 & 1.087093 & 1.087698 & 6.05E-04\\
4.0 & 1.501637 & 1.502164 & 5.26E-04 & 1.215055 & 1.214361 & 6.94E-04\\
4.5 & 1.974847 & 1.973759 & 1.09E-03 & 1.784686 & 1.785414 & 7.27E-04\\
5.0 & 1.705620 & 1.704236 & 1.38E-03 & 2.828001 & 2.827758 & 2.43E-04\\
\hline
\end{tabular}
}
\vspace*{0.2cm}
\caption{\it Numerical results for Example \ref{EX2} }\label{Table-ex2}
\end{table}
\begin{figure}[H]
\centering
\includegraphics[scale=.65,trim={.6cm 0.6cm 4cm 0.5cm},clip]{fig.eps}
\vspace*{0.1cm}
\caption{\it Comparison between   approximate solutions by the \emph{HBPLM} and \emph{RK4}}\label{Fig1}
\end{figure}
}
\end{example}
%*********** Example 3
\begin{example}\label{EX3}
\emph{In this example, the following system is considered as a linear stiff system of ODEs \cite{Biazar2015}:
\[%\begin{equation}\label{example3}
\left \{
\begin{aligned}
 u'_1(t)&=-u_1(t)+95\,u_2(t),\\
u'_2(t)&=-u_1(t)-97\,u_2(t).\\
\end{aligned}
\right.
\]%\end{equation}
For the initial conditions $u_1(0)=1$, $u_2(0)=1$, the exact solution is $u_1(t)=\frac{1}{47}(95e^{-2t}-48e^{-96t})$, $u_2(t)=\frac{1}{47}(48e^{-96t}-e^{-2t})$. The absolute errors of $u_1$ and $u_2$ by using the HBPLM for $N=10,\,M=12$, the HPM \cite{Biazar2015}, and the  rational HPM (RHPM) \cite{Biazar2015} are shown in Table \ref{Table-ex3}.
\begin{table}[h]
\centering
\resizebox{1\textwidth}{!}{
\begin{tabular}{lccc|ccc}
\hline
\multirow{2}{*}{t} & \multicolumn{3}{c}{Absolute error of $u_1$} &  \multicolumn{3}{c}{Absolute error of $u_2$}\\
\cline{2-4}
\cline{5-7}
& HBPLM & RHPM \cite{Biazar2015}& HPM \cite{Biazar2015}& HBPLM& RHPM  \cite{Biazar2015}& HPM  \cite{Biazar2015} \\
\hline
0.5 & 2.59E-16 & 1.98E-06 & 1.46E-06 & 5.07E-18 & 4.32E-05 &6.83E-05 \\
1.0 & 2.33E-16 & 8.71E-07& 3.73E-05 & 8.82E-17& 3.95E-05& 1.54E-03 \\
1.5 & 4.96E-18 & 1.22E-06&1.79E-04&1.36E-17&4.24E-05 &6.56E-03\\
2.0 & 1.55E-15&8.19E-07&4.53E-04&3.09E-17&3.11E-07 &1.45E-02\\
2.5 & 1.99E-15 & 1.54E-06& 8.13E-04& 4.82E-17& 2.94E-05& 2.25E-02\\
3.0 & 2.61E-17 & 7.77E-07& 1.18E-03& 5.52E-17& 1.36E-05& 2.75E-02\\
3.5 & 1.06E-16 & 4.38E-07& 1.48E-03& 1.27E-17& 8.66E-06& 2.79E-02\\
4.0 & 1.25E-16 & 1.39E-06& 1.68E-03& 2.03E-17& 1.82E-05& 2.41E-02\\
4.5 & 1.40E-16 & 1.87E-06& 1.75E-03& 2.54E-17& 1.65E-05& 1.73E-02\\
5.0 & 1.50E-16 & 1.93E-06& 1.73E-03& 2.98E-17& 1.02E-05& 9.20E-03\\
\hline
\end{tabular}
}
\vspace*{0.2cm}
\caption{\it Absolute errors for Example \ref{EX3} }\label{Table-ex3}
\end{table}
}
\end{example}
%********** example 4
\begin{example}\label{EX4}
\emph{As a non-linear stiff system of ODEs, we consider the following system \cite{Biazar2015}:
\[%\begin{equation}\label{example4}
\left \{
\begin{aligned}
 u'_1(t)&=-1002\,u_1(t)+1000\,u_2^2(t),\\
 u_2'(t)&=u_1(t)-u_2(t)-u_2^2(t),\\
\end{aligned}
\right.
\]%\end{equation}
with the initial conditions $u_1(0)=1$, $u_2(0)=1$, and the exact solution $u_1(t)=e^{-2t}$, $u_2(t)=e^{-t}$. The absolute errors of $u_1$ and $u_2$ using the HBPLM for $N=4,\,M=12$, the HPM \cite{Biazar2015} and the RHPM \cite{Biazar2015} are shown in Table \ref{Table-ex4}.
\begin{table}[h]
\centering
\resizebox{1\textwidth}{!}{
\begin{tabular}{lccc|ccc}
\hline
\multirow{2}{*}{t} & \multicolumn{3}{c}{Absolute error of $u_1$} &  \multicolumn{3}{c}{Absolute error of $u_2$}\\
\cline{2-4}
\cline{5-7}
& HBPLM & RHPM \cite{Biazar2015}& HPM \cite{Biazar2015}& HBPLM & RHPM \cite{Biazar2015}& HPM \cite{Biazar2015} \\
\hline
0.5 & 2.54E-12 & 3.00E-10 & 0.00E+00 & 1.09E-15 & 1.26E-08 & 0.00E+00 \\
1.0    & 6.30E-13 & 1.50E-08 & 0.00E+00 & 5.36E-16 & 1.49E-07 & 2.80E-09 \\
1.5 & 6.12E-14 & 2.75E-07 & 0.00E+00 & 5.52E-18 & 1.70E-07 & 1.75E-06\\
2.0    & 2.05E-13 & 4.61E-07 & 2.80E-09 & 3.11E-16 & 1.80E-07 & 1.66E-04\\
2.5 & 6.07E-13 & 3.69E-07 & 9.69E-08 & 7.05E-16 & 4.51E-08 & 5.62E-03\\
3.0    & 1.61E-14 & 6.80E-08 & 1.75E-06 & 3.87E-16 & 1.75E-07 & 9.93E-02\\
3.5 & 7.18E-15 & 3.99E-07 & 2.01E-05 & 2.83E-16 & 4.54E-08 & 1.12E+00\\
4.0    & 3.70E-15 & 8.09E-08 & 1.66E-04 & 7.72E-16 & 1.91E-07 & 9.09E+00\\
4.5 & 3.74E-15 & 5.19E-07 & 1.07E-03 & 3.88E-16 & 2.27E-07 & 5.75E+01\\
5.0    & 4.26E-14 & 1.19E-06 & 5.62E-03 & 2.12E-16 & 7.20E-07 & 2.98E+02\\
\hline
\end{tabular}
}
\vspace*{0.2cm}
\caption{\it Absolute errors for Example \ref{EX4} }\label{Table-ex4}
\end{table}}
\end{example}
%********** Example 5
\begin{example}\label{EX5}
\emph{Consider a system of ODEs of order two as \cite{Biazar2007}
 \begin{equation}\label{example5}
\left \{
\begin{aligned}
 u''_1(t)-u''_2(t)+u_1(t)-4u_2(t)&=0,\\
 u'_1(t)+u'_2(t)&=\cos t+2\cos 2t,\\
\end{aligned}
\right.
 \end{equation}
with the initial conditions $u_1(0)=0$, $u_2(0)=0$, $u'_1(0)=1$, $u'_2(0)=2$, and the exact solution $u_1(t)=\sin t$, $u_2(t)=\sin 2t$. We convert system \eqref{example5} into a system of ODEs of order one and solve it using the HBPLM. To do this, we define functions $w_1(t)$, $w_2(t)$, $w_3(t)$, and $w_4(t)$ as
\begin{equation*}
w_1(t)=u_1(t),\qquad w_2(t)=u'_1(t), \qquad w_3(t)=u_2(t),\qquad w_4(t)=u'_2(t).
\end{equation*}
Now, the following system is obtained:
\begin{equation}\label{convertedsysEx5}
\left \{
\begin{aligned}
 w'_1(t)&=w_2(t),\\
 w'_2(t)&=-\frac{1}{2}w_1(t)+2w_3(t)-\frac{1}{2}(\sin t+4\,\sin 2t),\\
 w'_3(t)&=w_4(t),\\
 w'_4(t)&=\frac{1}{2}w_1(t)-2w_3(t)-\frac{1}{2}(\sin t+4\,\sin 2t),\\
\end{aligned}
\right.
\end{equation}
where the initial conditions are $w_1(0)=0$, $w_2(0)=1$, $w_3(0)=0$, and $w_4(0)=2$. In Table \ref{Table-ex51}, the absolute errors of $u_1$ and $u_2$, using the HBPLM for $N=2,\,M=8$, and VIM \cite{Biazar2007} are compared. We point out that according to the definition of the functions $w_2$ and $w_4$, the approximate solutions of $u'_1$ and $u'_2$ are also obtained by solving   system \eqref{convertedsysEx5}. By considering $u'_1(t)=\cos t$, $u'_2(t)=2\cos 2t$, as the exact solutions of the derivatives, the absolute errors of $u'_1$ and $u'_2$ are shown in Table \ref{Table-ex52}. Compared to the errors obtained for $u_1$ and $u_2$ by the HBPLM, the errors for $u'_1$ and $u'_2$ are larger. The fact is that the numerical errors are usually larger for the derivatives.
\begin{table}[h]
        \begin{subtable}{.53\linewidth}
      \centering
        \caption{}\label{Table-ex51}
       \begin{tabular}{l|ccc|ccc|}
\hline
\multirow{2}{*}{t}& \multicolumn{3}{c}{$u_1$} &  \multicolumn{3}{c}{$u_2$}\\
\cline{2-4}
\cline{5-7}
& HBPLM & &VIM \cite{Biazar2007}&HBPLM & &VIM \cite{Biazar2007} \\
\hline
0.1 & 5.85E-10& & 4.10E-09  & 3.45E-10 && 7.40E-09\\
0.2 & 7.36E-10 && 9.60E-09  & 6.52E-10 && 6.60E-09\\
0.3 & 7.10E-10 && 2.50E-08  & 6.95E-10 && 1.20E-08\\
0.4 & 5.29E-10 && 2.50E-08  & 6.78E-10 && 2.30E-08\\
0.5 & 2.50E-10 && 2.70E-08  & 1.39E-10 && 2.90E-08\\
0.6& 6.25E-10 && 4.00E-08  & 1.17E-09 && 3.40E-08 \\
0.7 & 7.63E-10 && 3.50E-08  & 9.83E-10 && 2.80E-08\\
0.8 & 7.12E-10 && 4.50E-08  & 7.91E-10 && 3.90E-08\\
0.9& 5.03E-10 & &6.10E-08  & 2.67E-10 && 6.10E-08\\
1.0 & 4.12E-10 && 1.50E-08  & 4.15E-10 && 1.70E-07\\
\hline
\end{tabular}
    \end{subtable}
    \begin{subtable}{.6\linewidth}
      \centering
        \caption{}\label{Table-ex52}
        \begin{tabular}{|ccc}
\hline
\multicolumn{3}{c}{HBPLM} \\
\cline{1-3}
$u'_1$&&$u'_2$\\
\hline
1.01e-08&&1.10e-08\\
4.02e-09&&4.23e-09\\
6.48e-09&&6.59e-09\\
1.29e-08&&1.37e-08\\
3.89e-08&&4.04e-08\\
8.00e-09&&8.03e-09\\
2.56e-09&&2.62e-09\\
8.13e-09&&8.11e-09\\
1.52e-08&&1.52e-08\\
2.12e-08&&2.12e-08\\
\hline
\end{tabular}
    \end{subtable}
\vspace*{0.2cm}\caption{\it Absolute errors for \emph{\bf(a)} $u_1,\,u_2$ and \emph{\bf(b)} $u'_1,\,u'_2$ in Example \ref{EX5} }
\end{table}
}
\end{example}
%**********Example 6
\begin{example}\label{EX6}
\emph{As the last example, we consider Duffing's equation as \cite{Vahidi2018}
\begin{equation}\label{example6}
u''(t)+u'(t)+u(t)+u^3(t)=\cos^3 t-\sin t,\qquad u(0)=1,\, u'(0)=0.
\end{equation}
Considering $w_1(t)=u(t)$ and $w_2(t)=u'(t)$, system \eqref{example6} is converted into the following system of ODEs:
\begin{equation}\label{convertedsysEx6}
\left \{
\begin{aligned}
 w'_1(t)&=w_2(t),\\
w'_2(t)&=\cos^3 t-\sin t-w_1(t)-w_2(t)-w_1^3(t),
\end{aligned}
\right.
\end{equation}
with the initial conditions $w_1(0)=1$ and $w_2(0)=0$. The HBPLM is used for solving system \eqref{convertedsysEx6} for $N=2,\,M=8$. Then, the results obtained for the solutions of \eqref{convertedsysEx6} and consequently for equation \eqref{example6} with the exact solution $u(t)=\cos t$ are presented in Table \ref{Table-ex6}.
\begin{table}[h]
\centering
\resizebox{1\textwidth}{!}{
\begin{tabular}{lcccccc}
\hline
\multirow{1}{*}{t} &\multirow{1}{*}{\qquad\qquad\quad}&\multirow{1}{*}{Exact}&\multirow{1}{*}{\qquad\qquad\quad}& \multirow{1}{*}{HBPLM}&\multirow{1}{*}{\qquad\qquad\quad}& \multirow{1}{*}{Absolute error}\\
\hline
0.1 && 0.9950041652 && 0.9950041641 && 1.09E-09  \\
0.2 && 0.9800665778 && 0.9800665793 && 1.48E-09 \\
0.3 && 0.9553364891 && 0.9553364875 && 1.58E-09 \\
0.4 && 0.9210609940 && 0.9210609952 && 1.27E-09 \\
0.5 && 0.8775825618 && 0.8775825614 && 4.77E-10\\
0.6 && 0.8253356149 && 0.8253356156 && 7.00E-10 \\
0.7 && 0.7648421872 && 0.7648421864 && 8.50E-10\\
0.8 && 0.6967067093 && 0.6967067101 && 7.72E-10\\
0.9 && 0.6216099682 && 0.6216099677 && 5.30E-10 \\
1.0 && 0.5403023058 && 0.5403023063 && 4.51E-10 \\
\hline
\end{tabular}
}
\vspace*{0.2cm}\caption{\it Numerical results for Example \ref{EX6} }\label{Table-ex6}
\end{table}
}
\end{example}
\section{Conclusion}
This paper applied a hybrid functions method based on the combination of BPFs and Legendre polynomials to solve the systems of ODEs in linear and non-linear forms. By using the present method, a system of ODEs was reduced to a linear or non-linear system of algebraic equations. An upper bound of the error was obtained for the presented method. Various examples such as problem of HIV interactions with CD4$^+$ T cells, the predator and prey problem, and two examples of linear and non-linear stiff systems were solved. In addition, a system of ODEs of order two and  Duffing's equation were converted into the systems of ODEs of order one and then  solved. Simple application and the desired accuracy of the results obtained by this method show the efficiency of our proposed method.

%%%% Acknowledgment %%%%%%%%
\section*{Acknowledgement}
The authors are grateful to the editor and the anonymous referees for their valuable comments and suggestions, which have significantly improved this paper.
%%%% Bibliography  %%%%%%%%%%
\begin{thebibliography}{99}

%Type = Book
\bibitem{AbellBraselton}
{\sc M.\,L.\,Abell, J.\,P.\,Braselton}, {\em Introductory Differential Equations}, 5th ed., Academic Press, London, 2018.

%Type = Article
\bibitem{Asadyetal2005}
{\sc B.\,Asady, M.\,Tavassoli Kajani, A.\,Hadi Vencheh, A.\,Heydari}, {\em Solving second kind integral equations with hybrid Fourier and Block-pulse functions}, Appl. Math. Comput. {\bf160}(2005), 517--522.

%Type = Article
\bibitem{Biazar2015}
{\sc J.\,Biazar, M.\,A.\,Asadi, F.\,Salehi}, {\em Rational homotopy perturbation method for solving stiff systems of ordinary differential equations}, Appl. Math. Model. {\bf390}(2015), 1291--1299.

%Type = Article
\bibitem{Biazaretal2004}
{\sc J.\,Biazar, E.\,Babolian, R.\,Islam}, {\em Solution of the system of ordinary differential equations by Adomian decomposition method}, Appl. Math. Comput. {\bf147}(2004), 713--719.

%Type = Article
\bibitem{Biazar2007}
{\sc J.\,Biazar, H.\,Ghazvini}, {\em He's variational iteration method for solving linear and non-linear systems of ordinary differential equations}, Appl. Math. Comput. {\bf191}(2007), 287--297.

%Type = Article
\bibitem{BiazarMontazeri2005}
{\sc J.\,Biazar, R.\,Montazeri}, {\em A computational method for solution of the prey and predator problem}, Appl. Math. Comput. {\bf163}(2005) 841--847.

%Type = Book
\bibitem{Butcher}
{\sc J.\,C.\,Butcher}, {\em Numerical Methods for Ordinary Differential Equations}, 2nd ed., John Wiley \& Sons, 2008.

%Type = Book
\bibitem{canuto.etal}
{\sc C.\,Canuto, M.\,Y.\,Hussaini, A.\,Quarteroni, T.\,A.\,Zang}, {\em Spectral methods in fluid dynamics}, Springer-Verlag, Berlin, 1988.

%Type = Article
\bibitem{culshaw2000}
{\sc R.\,V.\,Culshaw, S.\,Ruan}, {\em A delay-differential equation model of HIV infection of CD4$^+$ T cells}, Math. Biosci. {\bf 165}(2000), 27--39.

%Type = Article
\bibitem{Darvishietal2007}
{\sc M.\,T.\,Darvishi, F.\,Khani, A.\,A.\,Soliman}, {\em The numerical simulation for stiff systems of ordinary differential equations}, Comput. Math. Appl. {\bf54}(2007), 1055--1063.

%Type = Article
\bibitem{Dubeyetal2001}
{\sc B.\, Dubey, B.\,Das, J.\,Hussain}, {\em A predator-prey interaction model with self and cross diffusion}, Ecol. Model. {\bf141}(2001), 67--76.

%Type = Article
\bibitem{essunger1994}
{\sc P.\,Essunger, A.\,S.\,Perleson}, {\em Modeling HIV infection of CD4$^+$ T-cell subpopulation}, J. Theoret. Biol. {\bf170}(1994), 367--391.

%Type = Article
\bibitem{FieldNoyes1974}
{\sc R.\,J.\,Field, R.\,M.\,Noyes}, {\em Oscillations in chemical systems. IV. Limit cycle behavior in a model of a real chemical reaction}, J. Chem. Phys. {\bf60}(1974), 1877--1884.

%Type = Article
\bibitem{Hesameddini2019}
{\sc E.\,Hesameddini, M.\,Riahi}, {\em Hybrid Legendre Block-pulse functions method for solving partial differential equations with non-local integral boundary conditions}, J. Inf. Optim. Sci. {\bf40}(2019), 1391--1403.

%Type = Article
\bibitem{HesameddiniShahbazi2018}
{\sc E.\,Hesameddini, M.\,Shahbazi}, {\em Hybrid Bernstein Block-pulse functions for solving system of fractional integro-differential equations}, Int. J. Comput. Math. {\bf95}(2018), 2287--2307.

%Type = Article
\bibitem{hsiao2009}
{\sc C.\,H.\,Hsiao}, {\em Hybrid function method for solving Fredholm and Volterra integral equations of the second kind}, J. Comput. Appl. Math. {\bf230}(2009), 59--68.

%Type = Article
\bibitem{kirschner1996}
{\sc D.\,Kirschner}, {\em Using mathematics to understand HIV immune dynamics}, Notices Amer. Math. Soc. {\bf43}(1996), 191--202.

%Type = Book
\bibitem{Lotka}
{\sc A.\,J.\,Lotka}, {\em Elements of Physical Biology}, William and Wilkins, Baltimore, 1925. (Reissued as Elements of  Mathematical Biology, Dover, New York, 1956)

%Type = Article
\bibitem{MaleknejadEbrahimzadeh2014}
{\sc K.\,Maleknejad, A.\,Ebrahimzadeh}, {\em An efficient hybrid pseudo-spectral method for solving
optimal control of Volterra integral systems}, Math. Commun. {\bf19}(2014), 417--435.

%Type = Article
\bibitem{MaleknejadHashemi2011}
{\sc K.\,Maleknejad, E.\,Hashemizadeh}, {\em Numerical solution of the dynamic model of a chemical reactor by hybrid functions}, Procedia. Comput. Sci. {\bf3}(2011), 908--912.

%Type = Article
\bibitem{MaleknejadSaeedipoor2017}
{\sc K.\,Maleknejad, E.\,Saeedipoor}, {\em An efficient method based on hybrid functions for Fredholm integral equation of the first kind with convergence analysis}, Appl. Math. Comput. {\bf304}(2017), 93--102.

%Type = Article
\bibitem{Marzbanetal2009}
{\sc H.\,R.\,Marzban, S.\,M.\,Hoseini, M.\,Razzaghi}, {\em Solution of Volterra’s population model via Block-pulse functions and Lagrange-interpolating polynomials}, Math. Methods Appl. Sci. {\bf 32}(2009), 127--134.

%Type = Article
\bibitem{MarzbanRazzaghi2006}
{\sc H.\,R.\,Marzban, M.\,Razzaghi}, {\em Solution of multi-delay systems using hybrid of Block-pulse functions and Taylor series}, J. Sound. Vib. {\bf292}(2006), 954--963.

%Type = Article
\bibitem{Marzbanetal2008}
{\sc H.\,R.\,Marzban, H.\,R.\,Tabrizidooz, M.\,Razzaghi}, {\em Hybrid functions for nonlinear initial-value problems with applications to Lane-Emden type equations}, Phys. Lett. A. {\bf372}(2008), 5883--5886.

%Type = Article
\bibitem{MashayekhiRazzaghi2016}
{\sc S.\,Mashayekhi, M.\,Razzaghi}, {\em Numerical solution of distributed order fractional differential equations by hybrid functions}, J. Comput. Phys. {\bf315}(2016), 169--181.

%Type = Article
\bibitem{merdan2007}
{\sc M.\,Merdan}, {\em Homotopy perturbation method for solving a model for HIV infection of CD4$^+$ T cells}, {\.I}stanbul Ticaret {\"U}niversitesi Fen Bilimleri Dergisi Yil, {\bf12}(2007), 39--52.

%Type = Article
\bibitem{merdan2011}
{\sc M.\,Merdan, A.\,G{\"o}kdo{\u g}an, A.\,Yildirim}, {\em On the numerical solution of the model for HIV infection of CD4$^+$ T cells}, Comput. Math. Appl. {\bf62}(2011), 118--123.

%Type = Article
\bibitem{Mirzaeeetal2019}
{\sc F.\,Mirzaee, S.\,Alipour, N.\,Samadyar}, {\em Numerical solution based on hybrid of Block-pulse and parabolic functions for solving a system of nonlinear stochastic Ito-Volterra integral equations of fractional order}, J. Comput. Appl. Math. {\bf349}(2019), 157--171.

%Type = Article
\bibitem{Mohammadietal2017}
{\sc F.\,Mohammadi, L.\,Moradi, D.\,Baleanu, A.\,Jajarmi}, {\em A hybrid functions numerical scheme for fractional optimal control problems: Application to nonanalytic dynamic systems}, J. Vib. Control. {\bf24}(2017), 5030--5043.

%Type = Article
\bibitem{nowak1996}
{\sc M.\,Nowak, C.\,R.\,M.\,Bangham}, {\em Population dynamics of immune responses to persistent viruses}, Science  {\bf272}(1996), 74--79.

%Type = Article
\bibitem{nowakl1991}
{\sc M.\,Nowak, R.\,May}, {\em Mathematical biology of HIV infections: Antigenic variation and diversity threshold}, Math. Biosci. {\bf106}(1991), 1--21.

%Type = Article
\bibitem{ongun2011}
{\sc M.\,Y.\,Ongun}, {\em The Laplace Adomian decomposition method for solving a model for HIV infection of CD4$^+$ T cells}, Math. Comput. Model. {\bf53}(2011), 597--603.

%Type = Article
\bibitem{perelson1990}
{\sc A.\,S.\,Perelson}, {\em Modeling the interaction of the immune system with HIV}, in:{\em Mathematical and Statistical Approaches to AIDS Epidemiology}, (C.\,Castillo-Chavez, Ed.), {\em Lecture Notes in Biomathematics}, Vol.\,83, Springer, Berlin, Heidelberg, 1990, 350--370.
%In: Castillo-Chavez C. (eds) Mathematical and Statistical Approaches to AIDS Epidemiology. Lecture Notes in Biomathematics, vol 83. Springer, Berlin, Heidelberg. https://doi.org/10.1007/978-3-642-93454-4_17


%Type = Article
\bibitem{perelson1993}
{\sc A.\,S.\,Perelson, D.\,E.\,Kirschner, R.\,D.\,Boer}, {\em Dynamics of HIV infection CD4$^+$ T cells}, Math. Biosci. {\bf114}(1993), 81--125.

%Type = Article
\bibitem{Rafieietal2017}
{\sc Z.\,Rafiei, B.\,Kafash, S.\,M.\,Karbasi}, {\em State-control parameterization method based on using hybrid functions of Block-pulse and Legendre polynomials for optimal control of linear time delay systems}, Appl. Math. Model. {\bf45}(2017), 1008--1019.

%Type = Article
\bibitem{Rahimietal2017}
{\sc J.\,Rahimi, D.\,D.\,Ganji, M.\,Khaki, K.\,Hosseinzadeh}, {\em Solution of the boundary layer flow of an Eyring-Powell non-Newtonian fluid over a linear stretching sheet by collocation method}, Alex. Eng. J. {\bf56}(2017), 621--627.

%Type = Article
\bibitem{Mahmoodetal2005}
{\sc A.\,Saad Mahmood, L.\,Casasus, W.\,Al-Hayani}, {\em The decomposition method for stiff systems of ordinary differential equations}, Appl. Math. Comput. {\bf167}(2005), 964--975.

%Type = Article
\bibitem{SahaSingh2018}
{\sc S.\,Saha Ray, S.\,Singh}, {\em Numerical solution of stochastic Volterra-Fredholm integral equations by hybrid Legendre Block-pulse functions}, Int. J. Nonlinear Sci. Numer. Simul. {\bf19}(2018), 1--9.

%Type = Article
\bibitem{Sahu2015}
{\sc P.\,K.\,Sahu, S.\,Saha Ray}, {\em Hybrid Legendre Block-pulse functions for the numerical solutions of system of nonlinear Fredholm-Hammerstein integral equations}, Appl. Math. Comput. {\bf270}(2015), 871--878.

%Type = Article
\bibitem{Singh2010}
{\sc V.\,K.\,Singh, R.\,K.\,Pandey, S.\,Singh}, {\em A stable algorithm for Hankel transforms using hybrid of Block-pulse and Legendre polynomials}, Comput. Phys. Commun. {\bf181}(2010), 1--10.

%Type = Article
\bibitem{srivastava2014}
{\sc V.\,K.\,Srivastava, M.\,K.\,Awasthi, S.\,Kumar}, {\em Numerical approximation for HIV infection of CD4$^+$ T cells mathematical model}, Ain Shams Engineering Journal {\bf5}(2014), 625--629.

%Type = Article
\bibitem{TavassoliVencheh2005}
{\sc M.\,Tavassoli Kajani, H.\,Vencheh}, {\em Solving second kind integral equations with hybrid Chebyshev and Block-pulse functions}, Appl. Math. Comput. {\bf163}(2005), 71--77.

%Type = Article
\bibitem{ThongmoonPusjuso2010}
{\sc M.\,Thongmoon, S.\,Pusjuso}, {\em The numerical solutions of differential transform method and the Laplace transform method for a system of differential equations}, Nonlinear Anal. Hybrid Syst. {\bf4}(2010), 425--431.

%Type = Article
\bibitem{Vahidi2018}
{\sc A.\,R.\,Vahidi, E.\,Babolian, Z.\,Azimzadeh}, {\em An improvement to the homotopy perturbation method for solving nonlinear Duffing’s equations}, Bull. Malays. Math. Sci. Soc. {\bf41}(2018), 1105--1117.

%Type = Article
\bibitem{lwang2006}
{\sc L.\,Wang, M.\,Y.\,Li}, {\em Mathematical analysis of the global dynamics of a model for HIV infection of CD4$^+$ T cells}, Math. Biosci. {\bf200}(2006), 44--57.

%Type = Article
\bibitem{yuzbasi2012}
{\sc S.\,Y\"{u}zba\c{s}{\i}}, {\em A numerical approach to solve the model for HIV infection of CD4$^+$ T cells}, Appl. Math. Model. {\bf36}(2012), 5876--5890.

%Type = Article
\bibitem{yuzbasi2016}
{\sc S.\,Y\"{u}zba\c{s}{\i}}, {\em An exponential collocation method for the solutions of the HIV infection model of CD4$^+$ T cells}, Int. J. Biomath. {\bf9}(2016), 1--15.

\end{thebibliography}
%\end{linenumbers}
\end{document} 