The high-order differential equations play a great important role in engineering science, such as mechanics, physics, static electricity, chemical reaction diffusion model, fluid dynamics and so on. They often occur in the form of multi-point boundary value problems (BVPs) for an $n$-th order ordinary differential equation. For example, an $m$-point BVP model of a dynamical system with $m$ degrees of freedom. Multi-point boundary value problems arise in a variety of physics area. Many problems in the theory of elastic stability can be handled by the multi-point problems (see [1]). Large size bridges are sometimes contrived with multi-point supports which correspond to a multi-point boundary value condition (see [2]). The existence and multiplicity of solutions of multi-point boundary value problems were studied by many authors, see [3-8] and the references therein. Dehghan and his research group proposed the Adomian decomposition method (ADM) (see [9]), variational iteration method (VIM) (see [10]), homotopy analysis method (HAM) (see [11]), sinc-collocation method (SCM) (see [12]) to solve the multi-point BVPs. Shooting method is used to solve multi-point BVPs in [13, 14]. Reproducing kernel method (RKM) is applied to solve the multi-point BVPs (see [15-19]). By combining the ADM and RKM, Geng and Cui proposed a method for solving nonlinear second order two-point BVP (see [20]). This method is also applied to solve fourth order three-point boundary value problem (see [21]). Differential transform method is used to deal with multi-point BVPs (see [22, 23]). Doha developed a method based on shifted Jacobi spectral method for high-order multi-point BVPs (see [24]). In [25], an efficient numerical algorithm based on Haar wavelet is presented to solve a class of linear and nonlinear nonlocal BVPs. To obtain a certain accuracy solution, it needs more collocation points. In [26], the authors applied the augmented block pulse function method for solving a system of arbitrary-order boundary value problem associated with initial conditions or multi-point boundary conditions in separated or non-separated forms. However, some of these methods are reliable and applicable for solving ordinary differential equations, most of them provide the solution only for a particular kind of differential equations or a particular kind of boundary conditions.

In this paper, we consider the multi-point boundary value problems in the following form

$ \mu^{(n)}(x)=g(x, \mu(x), \mu'(x), \cdots, \mu^{(n-1)}(x)), ~0\leq x\leq1 $

(1.1)

with the boundary conditions

$ \phi_{i}(\mu(x_{1}), \cdots, \mu(x_{m}), \mu'(x_{1}), \cdots, \mu'(x_{m}), \cdots, \mu^{(n-1)}(x_{1}), \cdots, \mu^{(n-1)}(x_{m}))=0, $

(1.2)

$i=1, 2, \cdots, n, $ where $x_{j}\in[0, 1]$, $j=1, 2, \cdots, m$. We assume that $g$ has the properties which guarantee the existence and uniqueness of the solution of problem (1.1) under conditions (1.2). The boundary conditions (1.2) can be more general form with linear and nonlinear cases.

Spectral methods play important roles in solving different kinds of differential equations. The main advantage of these methods lies in their accuracy for a given number unknowns. There are three widely used spectral methods, namely, Galerkin, collocation and Tau methods. During the last two decades, wavelets has been paid considerable attention from many scholars and has been applied in wide range of engineering disciplines. Especially, these wavelets constructed from orthogonal polynomials are widely used in seeking numerical solutions of various types of differential equations. Recently, the operational matrices for Legendre wavelet, Chebyshev wavelet and Bernoulli wavelet were extensively used to solve differential equations.

Motivated by the work mentioned above, the main aim of this paper is to find a simple and accurate method based on the second kind Chebyshev wavelets (SCW) for the solution of problem (1.1) under conditions (1.2). Also we want to further extend the applications of the second kind Chebyshev wavelets. Our method consists of reducing the multi-point BVPs to a set of algebraic equations through expanding the highest derivative $\mu^{(n)}(x)$ by the second kind Chebyshev wavelets with unknown coefficients. By solving these coefficients, we get the required approximate solution. In this paper, the application of the second kind Chebyshev wavelets spectral collocation method for finding an approximate solution for multi-point BVPs is investigated.

The rest part of this paper is organized as follows. In Section 2, the second kind Chebyshev wavelets and their properties are introduced. In Section 3, the proposed method is applied to solve high-order multi-point BVPs. In Section 4, we present some numerical examples to show the efficiency and applicability of this method for multi-point BVPs. Finally, a brief conclusion is given in Section 5.

The second kind Chebyshev wavelets $\psi_{n, m}(t)=\psi(k, n, m, t)$ have four arguments: $k$ can assume any positive integer, $n=1, 2, 3, \cdots, 2^{k-1}$, $m$ is the degree of the second kind Chebyshev polynomials. They are defined on the interval $[0, 1)$ as

$\begin{align} \psi_{n, m}(t)=\left\{ \begin{array}{ll} 2^{\frac{k}{2}}\widetilde{U}_{m}(2^{k}t-2n+1), & \frac{n-1}{2^{k-1}}\leq t < \frac{n}{2^{k-1}} \\ & \hbox{} \\ 0, & \hbox{otherwise, } \end{array} \right.\end{align} $

where

$ \widetilde{U}_{m}(t)=\sqrt{\frac{2}{\pi}}U_{m}(t), $

(2.1)

$m=0, 1, 2, \cdots, M-1$ and $M$ is a fixed positive integer. The coefficient in eq.(2.1) is for orthonormality, here $U_{m}(t)$ are the second kind Chebyshev polynomials of degree $m$ which are orthogonal with respect to the weight function $\omega(t)=\sqrt{1-t^{2}}$ on the interval $[-1, 1]$ and satisfy the following recursive formula

$ U_{0}(t)=1, \, \, U_{1}(t)=2t, U_{m+1}(t)=2tU_{m}(t)-U_{m-1}(t), \, \, m=1, 2, 3, \cdots. $

Note that when dealing with the second kind Chebyshev wavelets the weight function has to be dilated and translated as $\omega_{n}(t)=\omega(2^{k}t-2n+1).$

A function $f(x)\in L^{2}(\Bbb R)$ defined over $[0, 1)$ may be expanded by the second kind Chebyshev wavelets as

$ f(x)=\sum^{\infty}\limits_{n=1} \sum^{\infty}\limits_{m=0}c_{n, m}\psi_{n, m}(x), $

(2.2)

where $c_{n, m}=\langle f(x), \, \psi_{n, m}(x)\rangle_{L^{2}_{\omega}[0, 1)}=\displaystyle\int_{0}^{1}f(x)\psi_{n, m}(x)\omega_{n}(x)dx$ in which $\langle\cdot, \, \cdot\rangle_{L^{2}_{\omega}[0, 1)}$ denotes the inner product in $L^{2}_{\omega}[0, 1)$. If the infinite series in eq.(2.2) is truncated, then it can be written as

$ f(x)\cong\sum^{2^{k-1}}\limits_{n=1} \sum^{M-1}\limits_{m=0}c_{n, m}\psi_{n, m}(x)={\bf{C}}^{T}{\bf{\Psi}}(x), $

where ${\bf{C}}$ and ${\bf{\Psi}}(x)$ are $2^{k-1}M\times 1$ matrices given by

$ \begin{eqnarray}{\bf{C}}&=&\left( c_{1, 0} ~ c_{1, 1}~ \cdots~ c_{1, M-1}~c_{2, 0}~ c_{2, 1} ~ \cdots~ c_{2, M-1}~ \cdots ~c_{2^{k-1}, 0}~c_{2^{k-1}, 1}~ \cdots~ c_{2^{k-1}, M-1}\right)^{T}, \\ {\bf{\Psi}}(x)&=&\left( \psi_{1, 0}~\psi_{1, 1}~ \cdots~ \psi_{1, M-1}~ \psi_{2, 0}~ \psi_{2, 1}~ \cdots \psi_{2, M-1}~ \cdots~ \psi_{2^{k-1}, 0}~ \psi_{2^{k-1}, 1}~ \cdots~ \psi_{2^{k-1}, M-1} \right)^{T}.~~~~~~~~~ \end{eqnarray} $

(2.3)

The following theorem gives the convergence and accuracy estimation of the second kind Chebyshev wavelets expansion (see [28]).

Theorem 2.1  Let $f(x)$ be a second-order derivative square-integrable function defined on $[0, 1)$ with bounded second-order derivative, say $|f''(x)|\leq B$ for some constants $B$, then

(ⅰ) $f(x)$ can be expanded as an infinite sum of the second kind Chebyshev wavelets and the series converges to $f(x)$ uniformly, that is

$ f(x)=\sum\limits_{n=1}^{\infty}\sum\limits_{m=0}^{\infty}c_{n, m}\psi_{n, m}(x), $

where $c_{n, m}=\langle f(x), \, \psi_{n, m}(x)\rangle_{L^{2}_{\omega}[0, 1)}$.

(ⅱ)

$ \sigma_{f, k, M} < \frac{\sqrt{\pi}B}{2^{3}} \big(\sum\limits_{n=2^{k-1}+1}^{\infty}\frac{1}{n^{5}} \sum\limits_{m=M}^{\infty}\frac{1}{(m-1)^{4}}\big)^{\frac{1}{2}}, $

where

$ \sigma_{f, k, M}=\big(\int_{0}^{1}\big|f(x)-\sum\limits_{n=1}^{2^{k-1}} \sum\limits_{m=0}^{M-1}c_{n, m}\psi_{n, m}(x)\big|^{2}\omega_{n}(x)dx\big)^{\frac{1}{2}}. $

The operational matrix of integration of the second kind Chebyshev wavelets with omitted item has been derived in (see [28]), which plays very important role in solving high-order multi-point boundary value problems. For $k=3$ and $M=3$, let

$ {\bf{\Psi}}_{12}(t):=\left( \psi_{1, 0}~\psi_{1, 1} ~\psi_{1, 2} ~\psi_{2, 0} ~\psi_{2, 1}~\psi_{2, 2}~\psi_{3, 0}~ \psi_{3, 1} ~ \psi_{3, 2} ~ \psi_{4, 0}~ \psi_{4, 1} ~ \psi_{4, 1} \right)~^{T}, \\ {\bf{\widetilde{\Psi}}}_{12}(t):=\left( 0 ~0 ~ \psi_{1, 3} ~ 0 ~ 0~ \psi_{2, 3} ~ 0 ~ 0 ~ \psi_{3, 3} ~ 0~ 0 ~ \psi_{4, 3} \right)~^{T}. $

Then we have the following expression

$ \int_{0}^{t}{\bf{\Psi}}_{12}(s)ds={\bf{P}}_{12\times12}{\bf{\Psi}}_{12}(t)+{\bf{\widetilde{\Psi}}}_{12}(t), $

where

$\begin{align} {\bf{P}}_{12\times12}=\frac{1}{2^{3}}\left( \begin{array}{cccccccccccc} 1 & \frac{1}{2} & 0 & 2 & 0 & 0 & 2 & 0 & 0 & 2 & 0 & 0 \\ -\frac{3}{4} & 0 & \frac{1}{4} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \frac{1}{3} & -\frac{1}{6} & 0 & \frac{2}{3} & 0 & 0 & \frac{2}{3} & 0 & 0 & \frac{2}{3} & 0 & 0 \\ 0 & 0 & 0 & 1 & \frac{1}{2} & 0 & 2 & 0 & 0 & 2 & 0 & 0 \\ 0 & 0 & 0 & -\frac{3}{4} & 0 & \frac{1}{4} & 0& 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{3} & -\frac{1}{6} & 0 & \frac{2}{3} & 0 & 0 & \frac{2}{3} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & \frac{1}{2} & 0 & 2 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & -\frac{3}{4} & 0 & \frac{1}{4} &0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & \frac{1}{3} & -\frac{1}{6} & 0 & \frac{2}{3} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0& 0 & 0 & 0 &1 & \frac{1}{2} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0& 0 & 0 & 0 & -\frac{3}{4} & 0 & \frac{1}{4}\\ 0 & 0 & 0 & 0 & 0 & 0& 0 & 0 & 0 & \frac{1}{3} & -\frac{1}{6} & 0 \\ \end{array} \right).\end{align} $

In fact, the matrix ${\bf{P}}_{12\times12}$ can be written as

$\begin{align} {\bf{P}}_{12\times12}=\frac{1}{2^{3}}\left( \begin{array}{cccc} {\bf{L}}_{3\times3} & {\bf{F}}_{3\times3} & {\bf{F}}_{3\times3} & {\bf{F}}_{3\times3} \\ {\bf{0}}_{3\times3} & {\bf{L}}_{3\times3} & {\bf{F}}_{3\times3} & {\bf{F}}_{3\times3}\\ {\bf{0}}_{3\times3} & {\bf{0}}_{3\times3} & {\bf{L}}_{3\times3} & {\bf{F}}_{3\times3} \\ {\bf{0}}_{3\times3} & {\bf{0}}_{3\times3} & {\bf{0}}_{3\times3} & {\bf{L}}_{3\times3} \\ \end{array} \right), \end{align} $

where

$\begin{align} {\bf{L}}_{3\times3}=\left( \begin{array}{ccc} 1 & \frac{1}{2} & 0 \\ -\frac{3}{4} & 0 & \frac{1}{4} \\ \frac{1}{3} & -\frac{1}{6} & 0 \\ \end{array} \right), \, \, \, {\bf{F}}_{3\times3}=\left( \begin{array}{ccc} 2 & 0 & 0 \\ 0 & 0 & 0 \\ \frac{2}{3} & 0 & 0 \\ \end{array} \right) \, \, {\rm and}\, \, {\bf{0}}_{3\times3}=\left( \begin{array}{ccc} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0& 0 & 0 \\ \end{array} \right).\end{align} $

In general, when $M\geq4$, it follows

$ \int_{0}^{t}{\bf{\Psi}}(s)ds={\bf{P}}{\bf{\Psi}}(t)+{\bf{\widetilde{\Psi}}}(t), $

(2.4)

where ${\bf{\Psi}}(t)$ is given in (2.3) and ${\bf{P}}$ is a $2^{k-1}M\times2^{k-1}M$ matrix given by

$\begin{align} {\bf{P}}=\frac{1}{2^{k}}\left( \begin{array}{ccccc} {\bf{L}} & {\bf{F}} & {\bf{F}} & \cdots & {\bf{F}} \\ {\bf{0}} & {\bf{L}} & {\bf{F}} & \cdots & {\bf{F}} \\ \vdots & {\bf{0}} & \ddots & \ddots & \vdots \\ \vdots & \vdots & \vdots & \ddots & {\bf{F}} \\ {\bf{0}} & {\bf{0}} & \cdots& {\bf{0}} & {\bf{L}} \\ \end{array} \right), \end{align} $

here ${\bf{L}}$ and ${\bf{F}}$ are $M\times M$ matrices given by

$ \begin{eqnarray*}&&{\bf{L}}=\left( \begin{array}{ccccccc} 1 & \frac{1}{2} & 0 & 0 & \cdots & 0 & 0 \\ -\frac{3}{4} & 0 & \frac{1}{4} & 0 & \cdots & 0 & 0 \\ \frac{1}{3} & -\frac{1}{6} & 0 & \frac{1}{6} & \cdots & 0 & 0 \\ -\frac{1}{4} & 0 & -\frac{1}{8} & 0 & \ddots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \ddots & \ddots & \vdots \\ (-1)^{M-2}\frac{1}{M-1} & 0 & 0 & 0 & \ddots & 0 & \frac{1}{2(M-1)} \\ (-1)^{M-1}\frac{1}{M} & 0 & 0 & 0 & \cdots & -\frac{1}{2M} & 0 \\ \end{array} \right), \\ &&{\bf{F}}=\left( \begin{array}{ccccc} a_{1} & 0 & 0 & \cdots & 0 \\ a_{2} & 0 & 0 & \cdots & 0 \\ a_{3} & 0 & 0 & \cdots & 0 \\ \vdots & \vdots & \vdots & \cdots & \vdots \\ a_{M} & 0 & 0 & \cdots & 0 \\ \end{array} \right), \end{eqnarray*} $

where ${{a}_{i}}=\left\{ \begin{align} & \frac{2}{i},\ \ \ i\ \ \text{odd} \\ & \text{0,}\ \ \ \ i\ \ \text{even} \\ \end{align} \right.$ $i=1, 2, 3, \cdots, M$, and ${\bf{\widetilde{\Psi}}}(t)$ in eq.(2.4) is called as omitted item given by

$\begin{align} {\bf{\widetilde{\Psi}}}(t)=\frac{1}{2^{k}}\left( \begin{array}{ccccc} {\bf{L}}_{1} & {\bf{L}}_{2} &{\bf{L}}_{3}& \cdots & {\bf{L}}_{2^{k-1}} \\ \end{array} \right)^{T}, \end{align} $

where ${\bf{L}}_{i}$ are $1\times M$ matrices given by

$\begin{align} {\bf{L}}_{i}=\frac{1}{2M}\left( \begin{array}{cccccc} 0 & 0 & 0 & \cdots & 0 & \psi_{i, M} \\ \end{array} \right), \, \, \, i=1, 2, 3, \cdots, 2^{k-1}.\end{align} $

In this section, numerical solution of BVPs with multi-point boundary conditions based on the second kind Chebyshev wavelets is discussed. We assume that the highest derivative in the differential equation is approximated by the second kind Chebyshev wavelets as below

$ \mu^{(n)}(x)={\bf C}^{T}{\bf\Psi}(x), $

(3.1)

where ${\bf C}$ is an unknown vector which should be determined and ${\bf\Psi}(x)$ is the vector defined in (2.3). Before the further description of the proposed method, we introduce the following notations first

$ {\bf F}_{0}={\bf\widetilde{\Psi}}(x), \, \, \, {\bf F}_{n+1}(x)=\int_{0}^{x}{\bf F}_{n}(t)dt, \, \, n=0, 1, 2, \cdots,\\ {\bf F}_{n}(1)=\lim\limits_{x\rightarrow 1^{-}}{\bf F}_{n}(x), \, \, n=0, 1, 2, \cdots. $

Integrating eq.(3.1) and by eq.(2.4), we get

$ \begin{equation} \begin{aligned} \mu^{(n-1)}(x)=&{\bf C}^{T}\big({\bf P}{\bf\Psi}(x)+{\bf F}_{0}(x)\big)+\mu^{(n-1)}(0), \\ \mu^{(n-2)}(x)=&{\bf C}^{T}\big({\bf P}^{2}{\bf\Psi}(x)+{\bf P}{\bf F}_{0}(x)+{\bf F}_{1}(x)\big)+\mu^{(n-1)}(0)x+\mu^{(n-2)}(0), \\ \mu^{(n-3)}(x)=&{\bf C}^{T}\big({\bf P}^{3}{\bf\Psi}(x)+{\bf P}^{2}{\bf F}_{0}(x)+{\bf P}{\bf F}_{1}(x)+{\bf F}_{2}(x)\big)\\ &+\mu^{(n-1)}(0)\frac{x^{2}}{2}+\mu^{(n-2)}(0)x+\mu^{(n-3)}(0), \\ &\vdots\\ \mu(x)=&{\bf C}^{T}\big({\bf P}^{n}{\bf\Psi}(x)+\sum\limits_{i=0}^{n-1}{\bf P}^{n-i-1}{\bf F}_{i}(x)\big)+\sum\limits_{k=0}^{n-1}\frac{\mu^{(n-k-1)}(0)}{k!}x^{k}. \end{aligned} \end{equation} $

Without loss of generality, we may assume the values of $\mu^{(k)}(0)$, $k=0, 1, 2, \cdots, n-1$ are unknowns such that

$ \begin{equation} \begin{aligned} &\mu^{(n-1)}(x)={\bf C}^{T}\big({\bf P}{\bf\Psi}(x)+{\bf F}_{0}(x)\big)+d_{m}, \\ &\mu^{(n-2)}(x)={\bf C}^{T}\big({\bf P}^{2}{\bf\Psi}(x)+{\bf P}{\bf F}_{0}(x)+{\bf F}_{1}(x)\big)+d_{m}x+d_{m+1}(0), \\ &\mu^{(n-3)}(x)={\bf C}^{T}\big({\bf P}^{3}{\bf\Psi}(x)+{\bf P}^{2}{\bf F}_{0}(x)+{\bf P}{\bf F}_{1}(x)+{\bf F}_{2}(x)\big)+\frac{d_{m}}{2}x^{2}+d_{m+1}x+d_{m+2}(0), \\ &~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\vdots\\ &\mu(x)={\bf C}^{T}\big({\bf P}^{n}{\bf\Psi}(x)+\sum\limits_{i=0}^{n-1}{\bf P}^{n-i-1}{\bf F}_{i}(x)\big)+\sum\limits_{k=0}^{n-1}\frac{d_{m+k}}{k!}x^{k}, \end{aligned} \end{equation} $

where $d_{m+k}=\mu^{(n-k-1)}(0)$, $k=0, 1, 2, \cdots, n-1$. In order to calculate unknown vector ${\bf C}$ in (3.1), the following collocation points are considered

$ \xi_{l}=\frac{2l-1}{2^{k}M}, \, \, \, \, l=1, 2, 3, \cdots, 2^{k-1}M. $

We replace the expression of $\mu^{(i)}(x)$, $i=0, 1, 2, \cdots, n$, into systems of (1.1) and (1.2), and then substitute the collocation points as follows

$ \mu^{(n)}(\xi_{l})=g(\xi_{l}, \mu(\xi_{l}), \mu'(\xi_{l}), \cdots, \mu^{(n-1)}(\xi_{l})), ~~l=1, 2, 3, \cdots, 2^{k-1}M $

(3.2)

and

$ \phi_{i}\big(\mu(x_{1}), \cdots, \mu(x_{m}), \mu'(x_{1}), \cdots, \mu'(x_{m}), \cdots, \mu^{(n-1)}(x_{1}), \cdots, \mu^{(n-1)}(x_{m})\big)=0, \\ i=1, 2, \cdots, n. $

(3.3)

From eqs. (3.2) and (3.3), a system of $2^{k-1}M+n$ equations and $2^{k-1}M+n$ coefficients is obtained. Solving this system, we can get the unknown vector ${\bf C}$ and therefore the functions $\mu^{(i)}(x)$, $i=0, 1, 2, \cdots, n$, are identified. Note that we can do a few simple modifications when some of $\mu^{(i)}(0)$, $i=0, 1, 2, \cdots, n-1$ are given. Particularly if $\mu^{(i)}(0)$, $i=0, 1, 2, \cdots, n-1$ are all given, the system becomes an initial value problem and there is no need to consider any $d_{i}$, $i=m, m+1, \cdots, m+n-1$.

Remark 1  For linear multi-point BVPs, the coefficients $d_{i}$ are expressed by ${\bf\Psi}(x)$, ${\bf\widetilde{\Psi}}(x)$ and unknown vector ${\bf C}$. The unknown function $\mu(x)$ and its derivatives are substituted into eq.(3.2) and then collocation method is applied to generate linear system which can be solved using Gauss elimination method. The solution of this system gives us the values of unknown vector ${\bf C}$. We can calculate the approximate value of the unknown function $\mu(x)$ at any point using vector ${\bf C}$.

Remark 2  For nonlinear case, after substituting Chebyshev wavelet expression of $\mu(x)$ and its derivatives and collocation points into eqs. (3.2) and (3.3), a nonlinear system of $2^{k-1}M+n$ equations is obtained. This nonlinear system may be solved by using any iterative method for nonlinear system such as Newton's method, Broyden's method etc. Solving this system gives us the values of unknown vector ${\bf C}$ and $d_{i}$, which can be used to find approximate solution in a similar way as discussed for linear case.

Newton's iterative method for the nonlinear system ${\bf{F}}(x)=0$, where $ {\bf{F}}:\Omega \subset \mathcal{R}^n \rightarrow\mathcal{R}^n $ is defined by ${\bf{x}}_{k+1}={\bf{x}}_{k}-[{\bf{F}}^{'}({\bf{x}}_k)]^{-1} {\bf{F}}({\bf{x}}_k)$, where ${\bf{F}}^{'}({\bf{x}}_k)$ is the Jacobian matrix of point ${\bf{x}}_k$.

In this section, we give some numerical examples to demonstrate the efficiency and reliability of the proposed method. We also complete all the numerical computation and reported the absolute error at the selected points in tables and figures.

Example 1  Consider the following second-order four-point nonlinear ordinary differential equation [15, 29]

$ x(1-x)\mu''(x)+6\mu'(x)+2\mu(x)+\mu^{2}(x)=f(x), \, \, \, \, 0\leq x\leq1 $

with the boundary conditions

$ \mu(0)+\mu(\frac{2}{3})=\sinh(\frac{2}{3}), \, \, \mu(1)+\frac{1}{2}\mu(\frac{4}{5})=\frac{1}{2}\sinh(\frac{4}{5})+\sinh(1), $

where $f(x)=6\cosh(x)+\sinh(x)(2+x-x^{2}+\sinh(x))$. The exact solution is given by $\mu(x)=\sinh(x)$. In Table 1, we list the absolute errors at some different points obtained by the present method with $k=3$ and $M=4, 6, 8, 10$. As we see from this table, it is clear that the result obtained by the present method is superior to that obtained by the numerical methods given in [15, 29]. In Figure 1, we display the logarithmic plot for the maximal absolute errors at selected points with $k=3$ and $M=3$ through $10$ and the logarithmic plot for absolute errors with different $M$ and different methods.

Example 2  Consider the following third-order three-point nonlinear boundary value problems (see [13]) $\mu'''(x)=e^{-x}\mu^{2}(x), \, \, \, 0\leq x\leq1$ with the following nonlinear conditions

$ \begin{equation} \begin{aligned} &\mu(0)+2\mu'^{2}(\frac{1}{2})-\mu(1)-\mu''(1)-\mu(\frac{1}{2})=1-e^{\frac{1}{2}}, \\ &\mu'^{2}(0)-\mu''(0)\mu'(1)+\mu(\frac{1}{2})\mu''(\frac{1}{2})+\mu''(1)=1+e, \\ &\mu'(0)\mu(1)-\mu'(1)+\mu'(\frac{1}{2})+\mu(0)-\mu''(0)-\mu''(\frac{1}{2})=0. \end{aligned} \end{equation} $

The exact solution is $\mu(x)=e^{x}$. In Table 2, we compare the absolute errors at some different points obtained by the present method and the Newon-Broyden shooting method (NBSM) in [13]. In Figure 2, we show the logarithmic plot for the maximal absolute errors at selected points with $k=3$ and $M=3$ through $10$ and the logarithmic plot for absolute errors with different $M$ and different methods. It is evident from the figure that accuracy of the method increases with the increase of number $M$.

Example 3  Consider the following fourth-order three-point linear boundary value problem [18, 22]

$ \mu^{(4)}(x)=e^{x}\mu'''(x)-\mu(x)-e^{x}\cosh(x)+2\sinh(x)+1, \, \, \, 0\leq x\leq1 $

with boundary conditions given by

$ \mu(\frac{1}{4})=1+\sinh(\frac{1}{4}), \, \mu'(\frac{1}{4})=\cosh(\frac{1}{4}), \, \mu''(\frac{1}{4})=\sin(\frac{1}{4}), \\\, \mu(\frac{1}{2})-\mu(\frac{3}{4})=\sinh(\frac{1}{2})-\sinh(\frac{3}{4}). $

The exact solution is $\mu(x)=1+\sinh(x)$. In Table 3, we compare the absolute errors at some different points obtained by the present method and reproducing kernel methods in [18] and differential transform method in [22]. In Figure 3, we show the maximum absolute errors at selected points with $k=3$ and $M=3$ through $10$ and the logarithmic plot for absolute errors with different $M$ and different methods. Comparing with the Haar wavelet method in [25], we obtain more accurate solution in the case of the same number of collocation points. When the number of collocation points is 32, our result $L_{\infty}=5.8987\mbox{e-11}$, while in [25] $L_{\infty}=4.2510\mbox{e-7}$.

Example 4  Consider the following fifth-order four-point boundary value problems [18]

$ \mu^{(5)}(x)+\sin(2x)\mu'''(x)-\mu'(x)+\cos(2x)\mu(x)=-\sin x, \, \, \, \, 0\leq x\leq1 $

with boundary conditions given by

$ \begin{eqnarray*}&&\mu(0.1)=\sin(0.1), \, \, \mu'(0.1)=\cos(0.1), \, \, \mu(0.4)=\sin(0.4), \\ &&\mu'(0.4)=\cos(0.4), \, \, \mu(0.7)-\mu(0.9)=\sin(0.7)-\sin(0.9).\end{eqnarray*} $

The exact solution is $\mu(x)=\sin x$. In Table 4, we compare the absolute errors at some different points obtained by the present method and the reproducing kernel method (RKM) in [18]. In Figure 4, we show the maximum absolute errors at selected points with $k=3$ and $M=3$ through $10$ and the logarithmic plot for absolute errors with different $M$ and different methods.

Example 5  Consider following seventh-order three-point nonlinear boundary value problems [31]

$ \mu^{(7)}(x)=\mu(x)\mu'(x)-e^{x}(6+x-xe^{x}+x^{2}e^{x}), \, \, \, \, 0\leq x\leq1 $

with boundary conditions given by

$ \mu(0)=1, \, \, \mu(\frac{1}{2})=\frac{\sqrt{e}}{2}, \, \, \mu'(0)=0, \, \, \mu'(\frac{1}{2})=\\-\frac{\sqrt{e}}{2}, \, \, \mu''(0)=-1, \, \, \mu''(1)=-2e, \, \, \mu(1)=0. $

The exact solution is $\mu(x)=(1-x)e^{x}$. In Table 5, we compare the absolute errors at some different points obtained by the present method and variational iteration method (VIM) in [31]. In Figure 5, we show the maximum absolute errors at selected points with $k=3$ and $M=3$ through $10$ and the logarithmic plot for absolute errors with different $M$ and different methods.

Remark 3  In the all figures (left side), a logarithmic scale is used for the error-axis. It is obviously that the errors show an exponential decay, since in these semi-log representations one obverse that the error variations are essentially linear versus the polynomial degrees for all the examples. That is the so-called spectral accuracy as expected since the exact solution is smooth.

In this paper, the second kind Chebyshev wavelet collocation is applied to find numerical solutions of high-order multi-point boundary value problems. One of the advantages of the developed algorithm is that it can be applied on both linear and nonlinear problems and performs equally well in both cases. Another advantage is that high accurate approximate solutions can be achieved by using a few number of Chebyshev wavelet basis functions. Illustrative examples are presented to demonstrate the validity and applicability of the algorithm.