\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}
\AtBeginDocument{{\noindent\small
Ninth MSU-UAB Conference on Differential Equations and Computational
Simulations.
\emph{Electronic Journal of Differential Equations},
Conference 20 (2013), pp. 119--132.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2013 Texas State University - San Marcos.}
\vspace{9mm}}
\begin{document}\setcounter{page}{119}
\title[\hfilneg EJDE-2013/Conf/20/ \hfil Solving oscillatory problems]
{Solving oscillatory problems using a block hybrid
trigonometrically fitted method with two off-step points}
\author[F. F. Ngwane, S. N. Jator \hfil EJDE-2013/conf/20 \hfilneg]
{Fidele F. Ngwane, Samuel N. Jator} % in alphabetical order
\address{Fidele F. Ngwane \newline
Department of Mathematics,
University of south Caroline, Salkehatchie, SC 29488, USA}
\email{fifonge@yahoo.com}
\address{Samuel N. Jator \newline
Department of Mathematics and Statistics,
Austin Peay State University,
Clarksville, TN 37044, USA}
\email{Jators@apsu.edu}
\thanks{Published October 31, 2013.}
\subjclass[2000]{65L05, 65L06}
\keywords{First order system; hybrid method;
trigonometrically fitted method; \hfill\break\indent
off-step point; block method}
\begin{abstract}
A continuous hybrid method using trigonometric basis (CHMTB) with
two `off-step' points is developed and used to produce three discrete
hybrid methods which are simultaneously applied as numerical integrators
by assembling them into a block hybrid trigonometrically fitted method (BHTFM)
for solving oscillatory initial value problems (IVPs). The stability
property of the BHTFM is discussed and the performance of the method is
demonstrated on some numerical examples to show accuracy and efficiency
advantages.
\end{abstract}
\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{definition}[theorem]{Definition}
\newtheorem{example}[theorem]{Example}
\newtheorem{remark}[theorem]{Remark}
\allowdisplaybreaks
\section{Introduction}
In this article, we consider the first-order differential
equation
\begin{equation} \label{eq-1}
y' = f( x,y), \quad y(a) = y_{0}, \quad x \in [a,b],
\end{equation}
with oscillating solutions where $f:\mathbb{R} \times
\mathbb{R}^m\to \mathbb{R}^m$, $y, y_{0} \in \mathbb{R}^m$.
Oscillatory IVPs frequently arise in areas such as classical
mechanics, celestial mechanics, quantum mechanics, and biological
sciences. Several numerical methods based on the use of polynomial
basis functions have been developed for solving this class of
important problems (see Lambert \cite{Lam1,Lam2}, Hairer et al in
\cite{HWS}, Hairer \cite{HWS10}, and Sommeijer \cite{SJ}). Other methods based
on exponential fitting techniques which take advantage of the special
properties of the solution that may be known in advance have been
proposed (see Simos \cite{SIMOS}, Vanden et al \cite{VAN},
Vigo-Aguiar et al \cite{RAMOS}, Franco \cite{FRANCO},
Fang et al \cite{FANG}, Nguyen et al \cite{CONG}, and Jator et al \cite{JSF}).
The motivation governing the exponentially-fitted methods is inherent in
the fact that if the frequency or a reasonable estimate of it is
known in advance, these methods will be more advantageous than the
polynomial based methods.
The goal of this article is to construct a CHMTB which provides
three discrete methods that are combined and applied as a BHTFM which
takes the frequency of the solution as a priori knowledge. The
coefficients of the BHTFM are functions of the frequency and the
stepsize, hence the solutions provided by the proposed method are
exact if \eqref{eq-1} has periodic solutions with known frequencies. We are
motivated to use the hybrid method in order to increase the order of
the method, while preserving good stability properties. Hybrid
methods were first proposed to overcome the Dahlquist \cite{DQ}
barrier theorem whereby the conventional linear multistep method was
modified by incorporating off-step points in the derivation process
(see Gear \cite{GE1}, Gragg et al in \cite{GS}, Butcher
\cite{BU}, Gupta \cite{GUP}, Lambert \cite{Lam2}, and Kohfeld et al in
\cite{KT1}). These methods were shown to enjoy both higher
order and good stability properties, but included additional
off-grid functions. Gupta \cite{GUP} noted that the design of
algorithms for hybrid methods is more tedious due to the occurrence
of off-step functions which increase the number of predictors needed
to implement the methods. In order to avoid this deficiency, we
adopt a different approach based on a block-by-block implementation
instead of the traditional step-by-step implementation, generally
performed in a predictor-corrector mode.
Hence, we adopted the approach given in Jator et al in \cite{JSF},
where the CHMTB is used to generate a main and two additional
methods which are combined and used as a BHTFM to simultaneously
produce approximations
\[
\{y_{n+\mu}, y_{n+v}, y_{n+1}\}\text{ at a block of points }
\{x_{n+\mu}, x_{n+v}, x_{n+1}\},
\]
$h=x_{n+1}-x_n$, $n=0, \dots, N-1$, on a partition
$[a,b]$, where $\mu, v \in (0,1)$, $a, b \in {\mathbb R}$, $h$
is the constant stepsize, $n$ is a grid index and $N>0$ is the
number of steps. Block methods have also been considered by Shampine
and Watts \cite{SHW}. We emphasize that the BHTFM simultaneously
generates approximations $\{y_{n+\mu}, y_{n+v}, y_{n+1}\}$ to the exact
solutions $\{y(x_{n+\mu}), y(x_{n+v}), y(x_{n+1})\}$.
To apply the block method at the next block to obtain
$y_{n+2}$, the only necessary starting value is $y_{n+1}$, and the
loss of accuracy in $y_{n+1}$ does not affect subsequent points,
thus the order of the algorithm is maintained. It is unnecessary to
make a function evaluation at the initial part of the new block.
Thus, at all blocks except the first, the first function evaluation
is already available from the previous block.
The organization of this article is as follows.
In Section 2, we obtain a
trigonometric basis representation $U(x)$ for the exact solution
$y(x)$ which is used to generate two discrete methods which are
combined into a BHTFM for solving \eqref{eq-1}. The analysis and
implementation of the BHTFM are discussed in Section 3. Numerical
examples are given in Section 4 to show the accuracy and efficiency
of the BHTFM. Finally, we give some concluding remarks in Section 5.
\section{Development of method}
In this section, our objective is to construct a CHMTB which
produces three discrete methods as by-products. The main method has
the form
\begin{equation} \label{eq-2}
y_{n+1}=y_n + h (\beta_{0}(u) f_n + \beta_{1}(u) f_{n+1}
+ \beta_{v}(u) f_{n+v} + \beta_{\mu}(u) f_{n+\mu}),
\end{equation}
and the additional methods are given by
\begin{equation} \label{eq-3}
\begin{gathered}
y_{n+v}=y_n + h (\hat{\beta}_{0}(u) f_n + \hat{\beta}_{1}(u)f_{n+1}
+ \hat{\beta}_{v}(u) f_{n+v} + \hat{\beta}_{\mu}(u) f_{n+\mu}),
\\
y_{n+\mu}=y_n + h (\check{\beta}_{0}(u) f_n + \check{\beta}_{1}(u)f_{n+1}
+ \check{\beta}_{v}(u) f_{n+v} + \check{\beta}_{\mu}(u) f_{n+\mu}),
\end{gathered}
\end{equation}
where $u = w h$, $\beta_{j}(u)$, $\hat{\beta}_{j}(u)$, $\check{\beta}_{j}(u)$,
$\beta_{v}(u)$, $\hat{\beta}_{v}(u)$, $\check{\beta}_{v}(u)$, $\beta_{\mu}(u)$,
$\hat{\beta}_{\mu}(u)$ and $\check{\beta}_{\mu}(u)$, $j=0,1$, are
coefficients that depend on the stepsize and frequency. We note that
$y_{n+v}$ and $y_{n+\mu}$ are the numerical approximation to the analytical
solutions $y(x_{n+v})$, and $y(x_{n+\mu})$ respectively and
$f_{n+v}=f(x_{n+v},y_{n+v})$, $f_{n+\mu}=f(x_{n+\mu},y_{n+\mu})$,
$f_{n+j}=f(x_{n+j},y_{n+j})$ with $j = 0, 1$.
To obtain equations \eqref{eq-2} and \eqref{eq-3}, we
proceed by seeking to approximate the exact solution $y(x)$ on the
interval $[x_n, x_n+h]$ by the interpolating function $U(x)$ of
the form
\begin{equation} \label{eq-4}
U(x)= a_0 +a_1 x + a_2 x^2 + a_3 \sin(w x) + a_4 \cos(w x),
\end{equation}
where $a_0, a_1, a_2, a_3$ and $a_4$ are coefficients that must be
uniquely determined. We then impose that the interpolating function
equation \eqref{eq-4} coincides with the analytical solution at the
end point $x_n$ to obtain the equation
\begin{equation} \label{eq-5}
U(x_n)=y_n.
\end{equation}
We also require that the function \eqref{eq-4} satisfy
the differential equation \eqref{eq-1} at the points $x_{n+\mu}, x_{n+v},
x_{n+j}$, $j=0,1$ to obtain the following set of three equations:
%
\begin{equation} \label{eq-6}
U'(x_{n+\mu})=f_{n+\mu}, \quad U'(x_{n+v})=f_{n+v}, \quad
U'(x_{n+j})=f_{n+j}, \quad j=0,1.
\end{equation}
Equations \eqref{eq-5} and \eqref{eq-6} lead to a system of five
equations which is solved by Cramer's rule to obtain
$a_j$, $j = 0, 1, 2, 3, 4$. Our continuous CHMTB is constructed by
substituting the values of $a_j$ into equation \eqref{eq-4}.
After some algebraic
manipulation, the CHMTB is expressed in the form
\begin{equation} \label{eq-7}
U(x)=y_n + h(\beta_{0}(w, x) f_n + \beta_{1}(w, x) f_{n+1}
+ \beta_{v}(w, x) f_{n+v} + \beta_{\mu}(w, x) f_{n+\mu}),
\end{equation}
where $w$ is the frequency, $\beta_0(w,x)$,
$\beta_1(w,x)$, $\beta_v(w,x)$, and $\beta_\mu(w,x)$, are continuous
coefficients. The
continuous method \eqref{eq-7} is used to generate the main method
of the form \eqref{eq-2} and two additional methods of the form
\eqref{eq-3} by choosing, $v=\frac{1}{2}$ and $\mu=\frac{1}{4}$.
Thus, evaluating \eqref{eq-7} at
$x=\{x_{n+1}, x_{n+\frac{1}{2}}, x_{n+\frac{1}{4}}\}$ and letting
$u = w h$, we obtain the coefficients of \eqref{eq-2} and \eqref{eq-3} as
follows:
\begin{equation} \label{eq-8}
\begin{gathered}
\beta_0 = \frac{ \cos(\frac{u}{8}) (\csc(\frac{u}{4})^3)
\sin(\frac{u}{8}) (u - 2 \sin(\frac{u}{2})) }{2u} \\
\beta_1 = -( \frac{ \cos(\frac{u}{8}) \csc(\frac{u}{4})^3
\sin(\frac{u}{8}) (-u + 2 \sin(\frac{u}{2}))}{2u})\\
\beta_v = -( \frac{\cos(\frac{u}{8}) \csc(\frac{u}{4})^3
\sin(\frac{u}{8}) (u \cos(\frac{u}{2}) - 2 \sin(\frac{u}{2}))}{u}),
\end{gathered}
\end{equation}
and
\begin{equation} \label{eq-9}
\begin{gathered}
\hat{\beta}_0 = \frac{ (\csc(\frac{u}{8})^2) (u - 4 \sin(\frac{u}{4})) }{8 u}
\\
\hat{\beta}_v = \frac{(\csc(\frac{u}{8})^2) (u - 4 \sin(\frac{u}{4}))} {8 u}
\\
\hat{\beta}_\mu = -(\frac{ (\csc(\frac{u}{8})^2) (u \cos(\frac{u}{4})
- 4 \sin(\frac{u}{4}))}{4 u})
\\
\check{\beta}_0 = \frac{(\csc(\frac{u}{4})^3)
\sin(\frac{u}{8}) (8 u \cos(\frac{u}{8}) + 3 u \cos(\frac{(3 u)}{8})
- 8 (2 \sin(\frac{(3 u)}{8}) + \sin(\frac{(5 u)}{8})))} {16 u}
\\
\check{\beta}_1 = -(\frac{(\csc(\frac{u}{4})^3) (u \cos(\frac{u}{8})
- 8 \sin(\frac{u}{8})) \sin(\frac{u}{8})} {16 u})
\\
\check{\beta}_v = \frac{(3 + 3 \cos(\frac{u}{4}) + \cos(\frac{u}{2}))
( \csc(\frac{u}{4})^3) (u \cos(\frac{u}{8}) - 8 \sin(\frac{u}{8}))
\sin(\frac{u}{8})} {8 u}
\\
\check{\beta}_\mu = -( \frac{(\cos(\frac{u}{8})^2)
(\csc(\frac{u}{4})^3) \sin(\frac{u}{8}) (3 u \cos(\frac{u}{8})
+ 3 u \cos(\frac{(3 u)}{8}) - 16 \sin(\frac{(3 u)}{8}))} { 4 u}),
\end{gathered}
\end{equation}
with $\hat{\beta}_1 = 0 $ and $ \beta_\mu = 0$.
\section{Error analysis and stability}
\subsection{Local truncation error}
The Taylor series is used for small values of $u$ (see
Simos \cite{SIMOS}). Thus the coefficients in equation \eqref{eq-8}
can be expressed as
\begin{gather*} %10
\beta_0 = \frac{1}{6} + \frac{u^2}{720} + \frac{u^4}{80640}
+ \frac{u^6}{9676800} + \frac{u^8}{1226244096}
+ \frac{691 u^{10}}{111588212736000}+ \dots
\\
\beta_1 = \frac{1}{6} + \frac{u^2}{720} + \frac{u^4}{80640}
+ \frac{u^6}{9676800} + \frac{u^8}{1226244096}
+ \frac{691u^{10}}{111588212736000} + \dots
\\
\beta_v = \frac{2}{3} - \frac{u^2}{360} - \frac{u^4}{40320}
- \frac{u^6}{4838400} - \frac{u^8}{613122048} - \frac{691u^{10}}{55794106368000}
+ \dots
\end{gather*}
and the coefficients in equation \eqref{eq-9} can be expressed as
\begin{gather*} %11
\begin{aligned}
\hat{\beta}_0 &= \frac{1}{12} + \frac{u^2}{5760} + \frac{u^4}{2580480}
+ \frac{u^6}{1238630400} + \frac{u^8}{627836977152} \\
&\quad + \frac{691 u^{10}}{228532659683328000} + \dots
\end{aligned}\\
\begin{aligned}
\hat{\beta}_v &= \frac{1}{12} + \frac{u^2}{5760} + \frac{u^4}{2580480}
+ \frac{u^6}{1238630400} + \frac{u^8}{627836977152} \\
&\quad + \frac{691 u^{10}}{228532659683328000} + \dots
\end{aligned}\\
\begin{aligned}
\hat{\beta}_\mu &= \frac{1}{3} - \frac{u^2}{2880} - \frac{u^4}{1290240}
- \frac{u^6}{619315200} - \frac{u^8}{313918488576} \\
&\quad - \frac{691 u^{10}}{114266329841664000} + \dots
\end{aligned}\\
\begin{aligned}
\check{\beta}_0 &= \frac{37}{384} + \frac{67u^2}{184320}
+ \frac{401u^4}{165150720} + \frac{1649u^6}{79272345600}
+ \frac{1711u^8}{9132174213120} \\
&\quad + \frac{12094087 u^{10}}{7313045109866496000}
+ \dots
\end{aligned}\\
\begin{aligned}
\check{\beta}_1 &= \frac{1}{384} + \frac{13 u^2}{184320}
+ \frac{37 u^4}{33030144} + \frac{1091 u^6}{79272345600}
+ \frac{2087 u^8}{14350559477760} \\
&\quad + \frac{10191073 u^{10}}{7313045109866496000}
+ \dots
\end{aligned}\\
\begin{aligned}
\check{\beta}_v & = -\frac{7}{192} + \frac{7 u^2}{46080}
- \frac{11 u^4}{11796480} - \frac{29 u^6}{1415577600}
- \frac{12503 u^8}{50226958172160}\\
&\quad - \frac{659969 u^{10}}{261180182495232000}
+ \dots
\end{aligned}\\
\begin{aligned}
\check{\beta}_\mu & = \frac{3}{16} - \frac{3 u^2}{5120} - \frac{3 u^4}{1146880}
- \frac{31 u^6}{2202009600} - \frac{13 u^8}{155021475840} \\
&\quad - \frac{11747 u^{10}}{22571126882304000} + \dots.
\end{aligned}
\end{gather*}
Thus, the Local Truncation Errors (LTEs) for methods
\eqref{eq-2} and \eqref{eq-3} are given by
\begin{equation} %12
\begin{gathered}
\text{LTE}\eqref{eq-2} = -\frac{h^5}{2880}(w^2y^{(3)}(x_n) + y^{(5)}(x_n)),
\\
\text{LTE}\eqref{eq-3} = -\frac{h^5}{92160}(w^2y^{(3)}(x_n) + y^{(5)}(x_n)),
\\
\text{LTE}\eqref{eq-3} = -\frac{53h^5}{1474560}(w^2y^{(3)}(x_n) + y^{(5)}(x_n)).
\end{gathered}
\end{equation}
\begin{remark} \label{rmk3.1} \rm
We noticed that as $u \to 0$, the method (2)
reduces to the fourth-order method of Gragg and Stetter \cite{GS},
which is based on polynomial basis functions.
\end{remark}
\subsection{Stability}
The BHTFM is constructed by combining equations \eqref{eq-2} and
\eqref{eq-3}, where the coefficients of the method are explicitly
given by equations \eqref{eq-8} and \eqref{eq-9}. We then define the
block-by-block method for computing vectors $Y_0, Y_1, \cdots $ in
sequence (see \cite{FT}). Let the $\eta-$vector ($\eta = 3$ is the
number of points within the block) $Y_\gamma, Y_{\gamma -1}, F_\gamma,
F_{\gamma - 1}$ be given as
\begin{gather*}
Y_\gamma = (y_{n+ \frac{1}{4}}, y_{n+ \frac{1}{2}}, y_{n+1})^T, \quad
Y_{\gamma-1} = (y_{n - \frac{1}{4}}, y_{n - \frac{1}{2}}, y_n)^T, \\
F_\gamma = ( f_{n+ \frac{1}{4}}, f_{n+ \frac{1}{2}}, f_{n+1})^T, \quad
F_{\gamma-1} = (f_{n - \frac{1}{4}}, f_{n - \frac{1}{2}}, f_n)^T,
\end{gather*}
then the $1$-block $3$-point method for equation \eqref{eq-1} is
given by
\begin{equation}\label{eq-13}
Y_\gamma = \sum_{i=1}^1 A^{(i)}Y_{\gamma-i} + \sum_{i=0}^1
B^{(i)}F_{\gamma-i},
\end{equation}
where $A^{(i)}, B^{(i)}$, $i = 0, 1$ are $3 \times 3 $ matrices
whose entries are given by the coefficients of \eqref{eq-2} and \eqref{eq-3}.
\subsection*{Zero-stability}
\begin{definition}\label{df-1} \rm
The block method \eqref{eq-13} is zero stable provided the roots $R_j$,
$j=1, 2, 3$ of the first characteristic polynomial $\rho(R)$ specified
by
\begin{equation}\label{eq-14}
\rho(R) =\det \Big[ \sum_{i=0}^1 A^{(i)}R^{1-i} \Big] =
0, \quad A^{(0)} = -I,
\end{equation}
satisfies $|R_j| \leq 1$, $j= 1, 2, 3$ and for those roots with $|R_j| =1$,
the multiplicity does not exceed $1$ (see \cite{FT}).
\end{definition}
\subsection*{Consistency of BHTFM}
We note that the block method \eqref{eq-13} is consistent as it has
order $p>1$. We see from equation \eqref{eq-14} and definition \eqref{df-1}
that the block method \eqref{eq-13} is zero-stable since $\rho(R)=R^2(R-1)=0$
satisfies $|R_j| \leq 1$, $j=1, 2, 3$, and for those roots with
$|R_j| = 1$, the multiplicity does not exceed $1$. The block method
\eqref{eq-13} is thus convergent, as zero-stability plus consistency
equals convergence.
\subsection*{Linear stability of the BHTFM}
To analyze the linear stability of the BHTFM, we apply the method to the test
equation $y' =\lambda y$, where $\lambda$ is expected to run
through the (negative) eigenvalues of the Jacobian matrix. Then an
application of \eqref{eq-13} to the test equation yields
\begin{equation}\label{eq-15}
Y_\gamma = M(q;u)Y_{\gamma-i}, \quad q = h \lambda, \quad u = w h,
\end{equation}
where
\[
M(q;u) = (A^{(0)}- qB^{(0)})^{-1} (A^{(1)}+ qB^{(1)}),
\]
where the matrix $M(q;u)$ is the amplification matrix which determines
the stability of the method.
\begin{definition}\label{df-2} \rm
A region of stability is a region in the $q-u$ plane, throughout
which $|\rho(q;u)| \leq 1$, where $\rho(q;u)$ is the spectral radius of
$M(q;u)$.
\end{definition}
It is easily seen that the eigenvalues of $M(q;u)$ are
$\lambda_1 =0$, $\lambda_2 =0$, and
\begin{align*}
\lambda_3 &=
\Big(-16 (2 + q) (q^2 + u^2) \cos(\frac{u}{8}) + (6q^3 + 16u^2 + 6qu^2
+ q^2(16 + u^2)) \cos(\frac{3u}{8}) \\
&\quad + 16 q^2 \cos(\frac{5u}{8})
+ 10q^3 \cos(\frac{5u}{8}) + 16u^2 \cos(\frac{5u}{8})
+ 10qu^2 \cos(\frac{5u}{8}) \\
&\quad + 3q^2u^2\cos(\frac{5u}{8})
+ q^3u \sin(\frac{3u}{8}) + 3q^3u \sin(\frac{5 u}{8})\Big)\\
&\quad\div
\Big(16 (-2 + q) (q^2 + u^2) \cos(\frac{u}{8}) + (-10 q^3
+ 16 u^2 - 10 q u^2 \\
&\quad + q^2 (16 + 3 u^2)) \cos(\frac{3 u}{8}) + 16 q^2 \cos(\frac{5 u}{8})
- 6 q^3 \cos(\frac{5 u}{8}) + 16 u^2 \cos(\frac{5 u}{8}) \\
&\quad - 6 q u^2 \cos(\frac{5 u}{8})
+ q^2 u^2 \cos(\frac{5 u}{8})
- 3 q^3 u \sin(\frac{3 u}{8}) - q^3 u \sin(\frac{5 u}{8})\Big).
\end{align*}
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.4\textwidth]{fig1} % StaRegP3
\end{center}
\caption{The shaded region represents the truncated region of absolute
stability}
\label{StaReg}
\end{figure}
%
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.4\textwidth]{fig2} % ZeroPolesP3
\end{center}
\caption{$\lambda_{3}$ has three zeros$(\square)$ and no poles$(+)$ in $\mathbb{C}^-$, hence the BHTFM is $A_0$-stable.}
\label{ZerosPoles}
\end{figure}
We observed that in the $q-u$ plane the BHTFM is stable for
$ q \in [-100, 100]$, and $u \in [-\pi, \pi ]$.
Figure \ref{StaReg} is a plot of the stability region. Figure \ref{ZerosPoles}
shows the zeros and poles of $\lambda_3$.
\subsection{Implementation}
We emphasize that methods \eqref{eq-2} and \eqref{eq-3} are combined
to form the block method \eqref{eq-13}, which is implemented to solve
\eqref{eq-1} without requiring
starting values and predictors. For instance, if we let $n=0$ and
$\gamma = 1$ in \eqref{eq-13}, then
$(y_{\frac{1}{4}}, y_{\frac{1}{2}}, y_{1})^{T}$ are
simultaneously obtained over the sub-interval $[x_{0},x_{1}]$, as
$y_{0}$ is known from the IVP. Similarly, if $n = 1$ and $\gamma = 2$, then
$(y_{\frac{5}{4}}, y_{\frac{3}{2}}, y_{2} )^{T}$ are simultaneously obtained
over the sub-interval $[x_{1},x_{2}]$, as $y_{1}$ is known from the previous
block, and so on; until we reach the final sub-interval
$[x_{N-1},x_{N}]$. We note that for linear problems, we solve
\eqref{eq-1} directly using the feature solve[~~] in Matlab, while
nonlinear problems were solved by implementing the Newton's method
in Matlab enhanced by the feature fsolve[~~].
\section{Numerical examples}
In this section, we give numerical examples to illustrate the
accuracy (small errors) and efficiency (fewer number of function
evaluations (FNCs)) of the BHTFM. We find the approximate solution
on the partition $\pi_{N}$, where
\begin{equation}
\label{pi}
\pi_{N}:a=x_{0}