\documentclass[twoside]{article}
\usepackage{graphics}
\usepackage{graphicx}
\usepackage{latexsym}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{BSPMsample}
\usepackage{multirow,natlib,float}
%\usepackage {graphicx,natlib,fleqn,geometry,amsthm,amsmath,amsfonts,amssymb,float,bm,multirow}
\begin{document}

\title[ High Order Variable Mesh Size Exponential Finite Difference Method]
{High Order Variable Mesh Size Exponential Finite Difference Method for the Numerical Solutions of Two Point Boundary Value Problems }
%\thanks{The project is partially supported by...}}
%between [] goes the title that will appear on the top of every
%odd page, between {} goes the real title.
%\thanks{...} put your grant


\author[P. K. Pandey ]{P. K. Pandey }
%between [] goes the authors that will appear on the top of every
%even page, between {} goes the real author.

\address{P. K. Pandey\\Department of Mathematics, \\
Dyal Singh College (Univ. of Delhi), \\ Lodhi Road, New Delhi-110003, India\\
\\
\email:pramod\_10p@hotmail.com \\
}
\maketitle

\begin{abstract}
In this article, we presented a non-uniform mesh size high order exponential finite difference scheme for the numerical solutions of two point boundary value problems with  Dirichlet's boundary conditions. Under appropriate conditions, we have discussed the local truncation error and the convergence of the proposed method. Numerical experiments have been carried out to demonstrate the use and high order computational efficiency of the present method in several model problems. Numerical results showed that the proposed method is accurate and convergent. The order of accuracy is at least cubic which is in good agreement with the theoretically established order of the method.
\end{abstract}
%%put here Key words

\keywords  Exponential finite difference method, Cubic order, Singular Perturbed Problem, Two-point boundary value problem, High Order Non-uniform Mesh Size Method,  Nonlinear Method.


\tableofcontents


%% put here the subject class
\2000mathclass{65L05, 65L12}


\section{Introduction}\label{sec1}

In this article we considered a method for the numerical solution of the two-point boundary value problems of the form
\begin{equation}
  y''(x)=f(x,y),\quad a<x<b ,
\end{equation}
subject to the  boundary conditions
\begin{equation}
  y(a)=\alpha  \quad and \quad y(b)=\beta , \nonumber
\end{equation}
 where $\alpha $ and $\beta$  are real constants and $ f $ is continuous on $(x,y)$ for all $ x\in[a,b]$\quad $y \in\Re$. \\

Two point boundary value problems are of common occurrence in many areas of sciences and engineering. This class of problems has gained importance in the literature for the variety of their applications. In most cases it is impossible to obtain solutions of these problems using analytical methods which satisfy the given  boundary conditions. In these cases we resort to approximate solution of the problems and the last few decades have seen substantial progress in the development of approximate solutions of these problems. In all numerical solutions the continuous ordinary differential equation is replaced with a discrete approximation. In the literature, there are many different methods and approaches such as method of integration and discretization that are used to derive the approximate  solutions of these problems \cite{Ascher,Collatz,Keller,Stoer}. The finite difference method is one of several discretization techniques for obtaining approximate numerical solutions to problem (1).\\

The existence and uniqueness of the solution to problem (1) is assumed. We further assumed that problem (1) is well posed with continuous derivatives and that the solution depends differentially on the boundary conditions. The specific assumption on $f(x,y)$ to ensure existence and uniqueness will not be considered \cite{Keller,Stoer,Baxley}.\\

Over the last few decades, finite difference methods \cite{Roos,Hairer,Jain} have generated renewed interest and in recent years, variety of specialized techniques \cite{Van,Pandey-1} for the numerical solution of boundary value problems in ODEs have been reported in the literature. Recently, an exponential finite difference method with non-uniform step size was proposed in \cite{Pandey} for the numerical solution of two point boundary value problem. This method generated impressive numerical results for the problem (1). Hence, the purpose of this article is to propose an exponential finite difference method with non-uniform mesh size of higher order for problem (1). The main advantage for high order method is computationally more efficient than the second order method. In order to achieve a higher order method  more points are required for discretization approximation. But in proposed higher order method we have used only three points. \\

A method of at least cubic order is proposed for the numerical solution of linear boundary value problems (1). Our idea is to apply and extend the exponential finite difference method reported in \cite{Pandey-1} to discretize equation (1) in order to get a system of algebraic equations. In addition, if we apply a linearization technique and approximations, the method results in a tridiagonal matrix for the nodal values. The elements of this matrix depend on the source function i.e. right-hand side of the ordinary differential equation as well as on its partial derivatives with respect to the dependent variable and its first-order derivative.\\

We hope that others may find the proposed method an improvement and appealing to those existing finite difference methods for two-point boundary value problems.\\

We have presented our work in this article as follows. In the next section we derived a new variable mesh size exponential finite difference method. In Section 3, local truncation error and convergence of the new method are discussed in Section 4. The application of the proposed method to the problems in (1) has been presented and illustrative numerical results have been produced to show the efficiency of the new method in Section 5. Discussion and conclusion on the performance of the new method are presented in Section 6.

\section{The Variable Mesh Size Exponential Difference  Method }\label{sec2}
We define \,$N$ finite numbers of nodal points of the domain [\emph{a,b}], in which the solution of the problem (1) is desired, as $a\leq x_0<x_1<x_2<......<x_N<x_{N+1}=b,$ using nonuniform step length $\emph{h}_i$ such that  $ x_{i+1}=x_i+h_{i+1},\quad i=0,1,2,.....,N$ and $ r_i=\frac{h_{i+1}}{h_i}$. Suppose that we wish to determine the numerical approximation of the theoretical solution $y(x)$ of the problem (1) at the nodal point $x_i ,\quad i=1,2,.....,N$. We denote the numerical approximation of $y(x)$ at node $x=x_i$  as $y_i$ . Let us denote $f_i$ as the approximation of the theoretical value of the source function $f(x,y(x))$ at node $x=x_i ,\quad i=0,1,2,.....,N+1$. We can define other notations used in this article i.e. $f_{i\pm1}$,\quad and  $y_{i\pm1}$, in the similar way. Following the ideas in \cite{Pandey-1}, we propose an approximation to the theoretical solution $y(x_i)$ of the problem (1) by the exponential finite difference scheme as,
\begin{equation}
  a_2y_{i+1}+a_1y_i+a_0y_{i-1}=b_0h^2_if_i\exp(\phi(x_i)),\quad i=1,2,.......,N.
\end{equation}
where $a_0,a_1,a_2$ and $ b_0$ are unknown functions of the argument  $r_i$ and $\phi(x_i)$ is an unknown sufficiently differentiable function of \emph{x}.
Let us define a function $F_i(h,y)$ and associate it with (2) as,
\begin{equation}
  F_i(h,y)\equiv a_2y_{i+1}+a_1y_i+a_0y_{i-1}-b_0h^2_if_i\exp(\phi(x_i))=0 .
\end{equation}

Assume that $\phi(x_i)$ can be expanded in Taylor series about the point $x=x_{i-1}$. Hence we write $\phi(x_i)$ in Taylor series ,
\begin{equation}
  \phi(x_i)=\phi(x_{i-1})+h_i\phi'(x_{i-1})+\frac{h^2_i}{2}\phi''(x_{i-1})+ O(h^3_i).
\end{equation}
The application of (4) in the expansion of $\exp(\phi(x_i))$ will provide an $O(h^3_i)$ approximation of the form  as,
\begin{equation}
  \exp(\phi(x_i))=\exp(\phi(x_{i-1}))(1+h_i\phi'(x_{i-1})+\frac{h^2_i}{2}(\phi'^2(x_{i-1})+\phi''(x_{i-1})))+O(h^3_i)
\end{equation}
Expand  $F_i(h,y)$ in Taylor series about mesh point $x=x_i$ and using (5) in it, we have
\begin{multline}
 F_i(h,y)\equiv\{(a_0+a_1+a_2)y_i+h_i(r_ia_2-a_0)y'_i+\frac{h^2_i}{2}(r^2_ia_2+a_0)y''_i+\frac{h^3_i}{6}(r^3_ia_2-a_0)y^{(3)}_i+\\
 \frac{h^4_i}{24}(r^4_ia_2+a_0)y^{(4)}_i \} -b_0h^2_if_i\exp(\phi(x_{i-1}))(1+h_i\phi'(x_{i-1})+\frac{h^2_i}{2}(\phi'^2(x_{i-1})+\phi''(x_{i-1})))\\
 =0 .
\end{multline}
On comparing the coefficients of $ h^p_i ,p=0,1,2,3,4 $  both sides in (6), we  get the following system of nonlinear equations
\begin{equation}
  a_0+a_1+a_2=0,     \nonumber
\end{equation}
\begin{equation}
  r_ia_2-a_0=0,     \nonumber
\end{equation}
\begin{equation}
 ( r^2_ia_2+a_0)y''_i-2b_0f_i\exp(\phi(x_{i-1}))=0,     \nonumber
\end{equation}
\begin{equation}
  (r^3_ia_2-a_0)y^{(3)}_i-6b_0f_i\exp(\phi(x_{i-1}))\phi'(x_{i-1})=0 \nonumber
\end{equation}
\begin{equation}
  (r^4_ia_2+a_0)y^{(4)}_i-12b_0f_i\exp(\phi(x_{i-1}))(\phi'^2(x_{i-1})+\phi''(x_{i-1}))=0.
\end{equation}

To determine the unknown functions $a_0$, $a_1$, $b_0$ ,  $ \phi(x_{i-1}) $ , $ \phi'(x_{i-1})$  and  $ \phi''(x_{i-1})$  in (7), we have to assign arbitrary value to some unknown functions. To simplify the system of equations in (7), we have considered the following assumption:
\begin{equation}
 \phi(x_{i-1})=0 .
\end{equation}
Using (8) in (7) and solved the reduced system of equations, we obtained
\begin{equation}
  a_0=r_ia_2,                \nonumber
\end{equation}
\begin{equation}
  a_1=-(r_i+1)a_2,                \nonumber
\end{equation}
\begin{equation}
  b_0=\frac{r_i(r_i+1)a_2}{2},                \nonumber
\end{equation}
\begin{equation}
 \phi'(x_{i-1})=\frac{(r_i-1)y^{(3)}_i}{3f_i} \nonumber
\end{equation}
\begin{equation}
 \phi''(x_{i-1})=\frac{(r^2_i-r_i+1)y^{(4)}_i}{6f_i}-\frac{((r_i-1)y^{(3)}_i)^2}{(3f_i)^2}.
\end{equation}
Write $f'_i$ for $y^{(3)}_i$  and  $f''_i$ for $y^{(4)}_i$ in (9) and substituting the  values of $\phi(x_{i-1})$, $\phi'(x_{i-1} )$  and  $\phi''(x_{i-1} )$ from (8) and (9) in (4), we have
\begin{equation}
  \phi(x_i)=\frac{h_i(r_i-1)f'_i}{3f_i}+\frac{h^2_i}{2}\{ \frac{(r^2_i-r_i+1)f''_i}{6f_i}-\frac{((r_i-1)f'_i)^2}{(3f_i)^2}\}+O(h^3_i).
\end{equation}

Finally substitute the values of $a_0 $, $a_1$, $b_0$  and $\phi(x_i)$ from (9) and (10)  in  (2), we  obtain our proposed exponential finite difference method as
\begin{multline}
  y_{i+1}- (1+r_i)y_i+r_iy_{i-1}= \\
  \frac{h^2_i r_i(r_i+1)}{2}f_i\exp(\frac{h_i(r_i-1)f'_i}{3f_i}+\frac{h^2_i}{2}\{ \frac{(r^2_i-r_i+1)f''_i}{6f_i}-\frac{((r_i-1)f'_i)^2}{(3f_i)^2}\}).
\end{multline}

For each nodal point, we will obtain the nonlinear system of equations given by (11) or a linear system of equations if the source function is \emph{f(x)}. In the derived numerical method (11), the exponential function $\exp(\frac{h_i(r_i-1)f'_i}{3f_i}+\frac{h^2_i}{2}\{ \frac{(r^2_i-r_i+1)f''_i}{6f_i}-\frac{((r_i-1)f'_i)^2}{(3f_i)^2}\})$  has the argument $\frac{h_i(r_i-1)f'_i}{3f_i}+\frac{h^2_i}{2}\{ \frac{(r^2_i-r_i+1)f''_i}{6f_i}-\frac{((r_i-1)f'_i)^2}{(3f_i)^2}\}$. If $f_i$ in the denominator of the argument becomes zero in the domain of the solution, we  take the series expansion of the function $ \exp(\frac{h_i(r_i-1)f'_i}{3f_i}+\frac{h^2_i}{2}\{ \frac{(r^2_i-r_i+1)f''_i}{6f_i}-\frac{((r_i-1)f'_i)^2}{(3f_i)^2}\})$. Neglecting $O(h^3_i)$ and  higher order terms in the expansion, the method (11) becomes
\begin{equation}
  y_{i+1}- (1+r_i)y_i+r_iy_{i-1}=\frac{h^2_i}{24} r_i(r_1+1)(12f_i+4h_i(r_i-1)f'_i+h^2_i(r^2_i-r_i+1)f''_i).
\end{equation}
For the computational purpose in Section 4, let us define following $O(h^3_i)$ approximations :
\begin{equation}
  \overline{y}_{i+\frac{1}{2}}=c_{10}y_{i+1}+c_{11}y_i+c_{12}y_{i-1}+h^2_{i+1}c_{13}f_i ,
\end{equation}
where  $c_{10}=\frac{r^2_i-4}{8(r^2_i-1)}$, $c_{11}=\frac{3r^2_i+4r_i-4}{8(r_i-1)}$, $c_{12}=\frac{-3r^3_i}{8(r^2_i-1)}$ \quad and \quad $c_{13}=\frac{r_i+2}{16(r_i-1)}$.
\begin{equation}
  \overline{y}_{i-\frac{1}{2}}=c_{20}y_{i+1}+c_{21}y_i+c_{22}y_{i-1}+h^2_ic_{23}f_i ,
\end{equation}
where  $c_{20}=\frac{3}{8r_i(r^2_i-1)}$, $c_{21}=\frac{4r^2_i-4r_i-3}{8r_i(r_i-1)}$, $c_{22}=\frac{4r^2_i-1}{8(r^2_i-1)}$  $c_{23}=\frac{-(2r_i+1)}{16(r_i-1)}$.\\
Thus from (13)-(14), $ \overline{y}_{i\pm\frac{1}{2}} $ will provide $O(h^3_i)$ for  $ y_{i\pm\frac{1}{2}} $.\\
\\
Let us define following $O(h^5_i)$ approximations :
\begin{equation}
  \overline{\overline{y}}_{i+\frac{1}{2}}=c_{30}y_{i+1}+c_{31}y_i+c_{32}y_{i-1}+h^2_{i+1}(c_{33}f_{i+1}+c_{34}f_i) ,
\end{equation}
where  $c_{30}=\frac{3r^3_i-16r_i-8}{16(r^3_i-2r_i-1)}$, $c_{31}=\frac{5r^4_i+13r^3_i-16r_i-8}{16(r^3_i-2r_i-1)}$, $c_{32}=\frac{-5r^4_i}{16(r^3_i-2r_i-1)}$,\\   $c_{33}=\frac{r^3_i-7r_i-6}{96(1+2r_i-r^3_i)}$ \quad and \quad $c_{34}=\frac{4r^3_i+15r^2_i+17r_i+6}{96(r^3_i-2r_i-1)}$.
\begin{equation}
  \overline{\overline{y}}_{i-\frac{1}{2}}=c_{40}y_{i+1}+c_{41}y_i+c_{42}y_{i-1}+h^2_i(c_{43}f_i+c_{44}f_{i-1}) ,
\end{equation}
where  $c_{40}=\frac{5}{16r_i(r^3_i+2r^2_i-1)}$, $c_{41}=\frac{8r^4_i+16r^3_i-13r_i-5}{16r_i(r^3_i+2r^2_i-1)}$, $c_{42}=\frac{8r^4_i+16r^3_i-3r_i}{16r_i(r^3_i+2r^2_i-1)}$,\\
 $c_{43}=\frac{-(6r^3_i+17r^2_i+15r_i+4)}{96(r^3_i+2r^2_i-1)}$ \quad and \quad $c_{44}=\frac{-6r^3_i-7r^2_i+1}{96(r^3_i+2r_i-1)}$.\\
Thus from (15)-(16), $ \overline{\overline{y}}_{i\pm\frac{1}{2}} $ will provide $O(h^5_i)$ for  $ y_{i\pm\frac{1}{2}} $.\\
\\
Let define
\begin{equation}
  \overline{f}_{i+\frac{1}{2}}=f(x_{i+\frac{1}{2}},\overline{y}_{i+\frac{1}{2}}),\quad \overline{f}_{i-\frac{1}{2}}=f(x_{i-\frac{1}{2}},\overline{y}_{i-\frac{1}{2}})
\end{equation}
\begin{equation}
  \overline{\overline{f}}_{i+\frac{1}{2}}=f(x_{i+\frac{1}{2}},\overline{\overline{y}}_{i+\frac{1}{2}}),\quad and \quad  \overline{\overline{f}}_{i-\frac{1}{2}}=f(x_{i-\frac{1}{2}},\overline{\overline{y}}_{i-\frac{1}{2}}).
\end{equation}
Thus from (13),(14)  and (17), we find that $\overline{f}_{i\pm\frac{1}{2}}$ will provide $O(h^3_i)$ approximation for $f_{i\pm\frac{1}{2}}$. Similarly from (15), (16) and (18), we can find that $\overline{\overline{f}}_{i\pm\frac{1}{2}}$ will provide $O(h^5_i)$ approximation for $f_{i\pm\frac{1}{2}}$. \\

Using (17) and (18), let we define the following third order finite difference approximation for the terms  $ h_if'_i$ , $ h^2_if''_i$ in (11) and  (12) as :
\begin{equation}
  h_if'_i=c_{50}f_{i+1}+c_{51}f_{i-1}+c_{52}\overline{f}_{i+\frac{1}{2}}+c_{53}\overline{f}_{i-\frac{1}{2}},
\end{equation}
where $c_{50}=\frac{2-3r_i}{r_i(r_i+1)(2r_i+1)}$, $c_{51}=\frac{r_i(3-2r_i)}{(r_i+1)(r_i+2)}$, $c_{52}=\frac{4(3r_i-1)}{r_i(r_i+1)(r_i+2)}$  and \\ $c_{53}=\frac{4r_i(r_i-3)}{(r_i+1)(2r_i+1)}$.
\begin{equation}
  h^2_if''_i=c_{60}f_{i+1}+c_{61}f_i+c_{62}f_{i-1}+c_{63}\overline{\overline{f}}_{i+\frac{1}{2}}+c_{64}\overline{\overline{f}}_{i-\frac{1}{2}},
\end{equation}
where $c_{60}=\frac{4-6r_i}{r^2_i(r_i+1)(2r_i+1)}$, $c_{61}=\frac{2(4r^3_i-16r^2_i-5r_i+2)}{r^2_i(2r_i+1)}$, $c_{62}=\frac{2r_i(2r_i-3)}{(r_i+1)(r_i+2)}$, \\ $c_{63}=\frac{16(3r_i-1)}{r^2_i(r_i+1)(r_i+2)}$  and $c_{64}=\frac{16r_i(3-r_i)}{(r_i+1)(2r_i+1)}$.\\
Thus substituting the terms $h_if'_i$ and $h^2_if''_i$ from (19) and (20) in (11) or (12), we  obtained $O(h^3_i)$ method for computation of numerical solution of problems (1).

\section{Local Truncation Error }\label{sec4}
We can write the following expression for the term in (11) :
\begin{multline}
\exp(\frac{h_i(r_i-1)f'_i}{3f_i}+\frac{h^2_i}{2}\{ \frac{(r^2_i-r_i+1)f''_i}{6f_i}-\frac{((r_i-1)f'_i)^2}{(3f_i)^2}\})=1+\frac{h_i(r_i-1)f'_i}{3f_i} \\ +\frac{h^2_i(r^2_i-r_i+1)f''_i}{12f_i}+\frac{h^3_if'_i}{324f^3_i}(9(r^3_i+1)f_if''_i-4(r_i-1)^3(f'_i)^2)+O(h^4_i).
\end{multline}
From (11) and (21), the truncation error $T_i$ at the nodal point $ x=x_i$   may be written as \cite{ Jain,Varga,Golub},
\begin{multline}
  T_i=y_{i+1}-(1+r_i)y_i+r_iy_{i-1}-\frac{h^2_i(r^2_i+r_i)f_i}{2}(1+\frac{h_i(r_i-1)f'_i}{3f_i}+ \\
  \frac{h^2_i(r^2_i-r_i+1)f''_i}{12f_i}+\frac{h^3_if'_i}{324f^3_i}(9(r^3_i+1)f_if''_i-4(r_i-1)^3(f'_i)^2))+O(h^6_i). \nonumber
\end{multline}
By the Taylor series expansion of $y$ at nodal point $x=x_i$ and using $y''_i=f_i$ , $ y^{(3)}_i=f'_i $ and $ y^{(4)}=f''_i$  etc., we have
\begin{multline}
 T_i =(\frac{h^5_{i+1}}{120}-\frac{r_ih^5_i}{120})y^{(5)}_i+\frac{h^5_ir_i(r_i+1)y^{(3)}_i}{648(y''_i)^3}(9(r^3_i+1)y''_iy^{(4)}_i-4(r_i-1)^3(y^{(3)}_i)^2))+O(h^6_i)\\
 =\frac{h^5_i}{3240}(27r_i(r^4_i-1)y^{(5)}+\frac{5r_i(r_i+1)y^{(3)}_i}{(y''_i)^3}(9(r^3_i+1)y''_iy^{(4)}_i-4(r_i-1)^3(y^{(3)}_i)^2))+O(h^6_i).
\end{multline}
Thus from (22), $T_i$ can be written as :
\begin{equation}
 T_i =O(h^5_i).
\end{equation}
Thus we have obtained a truncation error at each node of  $O(h^5_i)$.

\section{Convergence of the Method}\label{sec5}
Let us substitute (19) and (20) into (12) and then simplify , we have
\begin{multline}
 y_{i+1}-(1+r_i)y_i+r_iy_{i-1}=\frac{h^2_ir_i(r_i+1)}{24}((12+(r^2_i-r_i+1)c_{61})f_i+\\
 (4(r_i-1)c_{50}+(r^2_i-r_i+1)c_{60})f_{i+1}+(4(r_i-1)c_{51}+(r^2_i-r_i+1)c_{62})f_{i-1}\\
 +4(r_i-1)c_{52}\overline{f}_{i+\frac{1}{2}}+(r^2_i-r_i+1)c_{63}\overline{\overline{f}}_{i+\frac{1}{2}}+4(r_i-1)c_{53}\overline{f}_{i-\frac{1}{2}}+(r^2_i-r_i+1)c_{64}\overline{\overline{f}}_{i-\frac{1}{2}}).  \nonumber
\end{multline}
Thus
\begin{multline}
 -y_{i+1}+(1+r_i)y_i-r_iy_{i-1}+\frac{h^2_i}{24}(\alpha_if_i+\gamma_if_{i+1}+\beta_if_{i-1}+\overline{\delta}_i\overline{f}_{i+\frac{1}{2}}+\\
 \overline{\overline{\delta}}_i\overline{\overline{f}}_{i+\frac{1}{2}}+\overline{\lambda}_i\overline{f}_{i-\frac{1}{2}}+\overline{\overline{\lambda}}_i\overline{\overline{f}}_{i-\frac{1}{2}})=0,
\end{multline}
where $\alpha_i=r_i(r_i+1)(12+(r^2_i-r_i+1)c_{61})$,\quad \\
$\beta_i=r_i(r_i+1)(4(r_i-1)c_{50}+(r^2_i-r_i+1)c_{60})$,\quad \\
$\gamma_i=r_i(r_i+1)(4(r_i-1)c_{51}+(r^2_i-r_i+1)c_{62})$,\quad $\overline{\delta}_i=4r_i(r^2_i-1)c_{52}$,\\
 $\overline{\overline{\delta}}_i=(r^2_i-r_i+1)c_{63}$,\quad $\overline{\lambda}_i=4r_i(r^2_i-1)c_{53}$,\quad and
 $\overline{\overline{\lambda}}_i=(r^2_i-r_i+1)c_{64}$.\\
Let us define
\begin{multline}
  \phi_1=\frac{h^2_1}{24}(\alpha_1f(x_1,y_1)+\gamma_1f(x_2,y_2)+\overline{\delta}_1f(x_{\frac{3}{2}},\overline{y}_{\frac{3}{2}})+\overline{\overline{\delta}}_1f(x_{\frac{3}{2}},\overline{\overline{y}}_{\frac{3}{2}})\\
 +\overline{\lambda}_1f(x_{\frac{1}{2}},\overline{y}_{\frac{1}{2}})+\overline{\overline{\lambda}}_1f(x_{\frac{1}{2}},\overline{\overline{y}}_{\frac{1}{2}}))+ \frac{h^2_1}{24}\beta_1f(x_0,y_0)+r_1y_0,\quad i=1\nonumber
\end{multline}
\begin{multline}
  \phi_i=\frac{h^2_i}{24}(\alpha_if(x_i,y_i)+\gamma_if(x_{i+1},y_{i+1})+\beta_if(x_{i-1},y_{i-1})+\overline{\delta}_if(x_{i+\frac{1}{2}},\overline{y}_{i+\frac{1}{2}})\\
 +\overline{\overline{\delta}}_if(x_{i+\frac{1}{2}},\overline{\overline{y}}_{i+\frac{1}{2}})
 +\overline{\lambda}_if(x_{i-\frac{1}{2}},\overline{y}_{i-\frac{1}{2}})+\overline{\overline{\lambda}}_if(x_{i-\frac{1}{2}},\overline{\overline{y}}_{i-\frac{1}{2}})),
  \quad 2\leq i \leq N-1 \nonumber
\end{multline}
\begin{multline}
  \phi_N=\frac{h^2_N}{24}(\alpha_Nf(x_N,y_N)+\beta_Nf(x_{N-1},y_{N-1})+\overline{\delta}_Nf(x_{N+\frac{1}{2}},\overline{y}_{N+\frac{1}{2}})\\
 +\overline{\overline{\delta}}_Nf(x_{N+\frac{1}{2}},\overline{\overline{y}}_{N+\frac{1}{2}})
 +\overline{\lambda}_Nf(x_{N-\frac{1}{2}},\overline{y}_{N-\frac{1}{2}})+\overline{\overline{\lambda}}_Nf(x_{N+\frac{1}{2}},\overline{\overline{y}}_{N+\frac{1}{2}}))\\
 +\frac{h^2_N}{24}\gamma_N f(x_{N+1},y_{N+1})+y_{N+1}.\quad i=N \nonumber
\end{multline}
Let us define column matrix $\boldsymbol{\phi}_{N\times1}$ and $\textbf{y}_{N\times1}$as
\begin{equation}
  \boldsymbol{\phi}=[\phi_1,\phi_2,............,\phi_N]'_{1\times N},\quad  \textbf{y}= [y_1,y_2,............,y_N]'_{1\times N},    \nonumber
\end{equation}
where $[.....]'$ is the transpose of a column matrix.

The difference method (24) represents a system of nonlinear equations in unknown $y_i,i=1,2,...,N$. Let us  write (24) in matrix form as,
\begin{equation}
  \textbf{Dy}+\boldsymbol{\phi}(\textbf{y}) = \boldsymbol{0},
\end{equation}
where \begin{equation}
  \textbf{D}=\left(
               \begin{array}{ccccc}
                  1+r_1 & -1 &  &  & 0 \\
                 -r_2 & 1+r_2 & -1 &  &  \\
                  &   -r_3 & 1+r_3 & -1 &  \\
                 .. & .. & .. & .. & .. \\
                 .. & .. & .. & .. & .. \\
                 0 &  &  & -r_N & 1+r_N \\
               \end{array}
             \right)_{N\times N}                  \nonumber
\end{equation}
is a tridiagonal matrix. Let \textbf{Y} be the exact solution of (24), so it will satisfy the matrix equation
\begin{equation}
  \textbf{DY}+\boldsymbol{\phi}(\textbf{Y})+\textbf{T} = \boldsymbol{0},
\end{equation}
where $\textbf{Y}$ is a column matrix of order $N\times1$ which can be obtained by replacing $y$ with $Y$ in matrix $\textbf{y}$  and $\textbf{T}$ is a truncation error matrix in which each element has $ O(h^5_i)$.
Let us define
\begin{equation}
  F_{i+\frac{1}{2}}=f(x_{i+\frac{1}{2}},Y_{i+\frac{1}{2}}) ,\quad F_{i-\frac{1}{2}}=f(x_{i-\frac{1}{2}},Y_{i-\frac{1}{2}}) ,            \nonumber
\end{equation}
\begin{equation}
  F_{i+1}=f(x_{i+1},Y_{i+1}) ,\quad f_{i+1}=f(x_{i+1},y_{i+1}),\quad F_{i-1}=f(x_{i-1},Y_{i-1}) ,            \nonumber
\end{equation}
\begin{equation}
  f_{i-1}=f(x_{i-1},y_{i-1}),\quad F_i=f(x_i,Y_i) , and \quad  f_i=f(x_i,y_i).     \nonumber
\end{equation}
After linearization of $\overline{f}_{i+\frac{1}{2}}$ and using (13) into it, we have
\begin{equation}
 \overline{f}_{i+\frac{1}{2}}= F_{i+\frac{1}{2}}+(\overline{y}_{i+\frac{1}{2}}-Y_{i+\frac{1}{2}})G_{i+\frac{1}{2}}
\end{equation}
where $G_{i+\frac{1}{2}}=(\frac{\partial f}{\partial Y})_{i+\frac{1}{2}}$ . But $\overline{y}_{i+\frac{1}{2}}$ provide $O(h^3_i)$ approximation for $Y_{i+\frac{1}{2}}$. Thus from (27), we have
\begin{equation}
  \overline{f}_{i+\frac{1}{2}}-F_{i+\frac{1}{2}} = O(h^3_i).
\end{equation}
Similarly, we can linearize $\overline{\overline{f}}_{i+\frac{1}{2}}$, $\overline{f}_{i-\frac{1}{2}}$, $\overline{\overline{f}}_{i-\frac{1}{2}}$, $f_{i+1}$, $f_{i-1}$, $f_i$ and use (14-16),  to obtain the following results :
\begin{equation}
  \overline{\overline{f}}_{i+\frac{1}{2}}-F_{i+\frac{1}{2}} = O(h^5_i).
\end{equation}
\begin{equation}
  \overline{f}_{i-\frac{1}{2}}-F_{i-\frac{1}{2}} = O(h^3_i).
\end{equation}
\begin{equation}
  \overline{\overline{f}}_{i-\frac{1}{2}}-F_{i-\frac{1}{2}} = O(h^5_i).
\end{equation}
\begin{equation}
  f_{i+1}-F_{i+1} = (y_{i+1}-Y_{i+1})G_{i+1}.
\end{equation}
\begin{equation}
 f_{i-1}-F_{i-1} = (y_{i-1}-Y_{i-1})G_{i-1},
\end{equation}
\begin{equation}
  f_i-F_i = (y_i-Y_i)G_i .
\end{equation}
From the difference of (25) and (26),  by taking the Taylor series expansion of $G_{i\pm1}$  about $x=x_i$. Neglecting the terms $O(h^3_i)$ and higher orders, we can write
\begin{equation}
  \boldsymbol{\phi}(\textbf{y})-\boldsymbol{\phi}(\textbf{Y})=\textbf{PE},
\end{equation}
where $\textbf{P}=(P_{lm})_{N\times N}$ is a tri-diagonal matrix  defined as
\begin{equation}
  P_{lm}=\frac{h^2_i}{24}(\alpha_iG_i) ,\quad     i=l=m ,\quad l=1,2,...,N ,   \qquad\qquad  \nonumber
\end{equation}
\begin{equation}
  P_{lm}=\frac{h^2_i}{24}\gamma_i(G_i+h_{i+1}(\frac{\partial G}{\partial x})_i), \quad m=l+1 ,\quad i=l=1,2,....,N-1,     \nonumber
\end{equation}
\begin{equation}
  P_{lm}=\frac{h^2_i}{24}\beta_i(G_i-h_i(\frac{\partial G}{\partial x})_i), \quad i=l=m+1 ,\quad m=1,2,....,N-2,   \nonumber
\end{equation}
and $\textbf{E}=[E_1,E_2,........,E_N]'_{ 1\times N }$ , where  $ E_i=(y_i-Y_i),i=1,2,....,N $ .

Let us assume that the solution of difference equation (11) has no roundoff error. So from (25), (26) and (35) we have
\begin{equation}
  (\textbf{D}+\textbf{P})\textbf{E}=\textbf{JE}=\textbf{T}.
\end{equation}
Let us define $G^0=\{G_i : i=1,2,...,N\}$,
\begin{equation}
 G_*= \min_{\substack {x\in[a,b]}}\frac{\partial f}{\partial Y}, and \quad G^*= \max_{\substack {x\in[a,b]}}\frac{\partial f}{\partial Y}, \nonumber
\end{equation}
such that
\begin{equation}
  0\leq G_*\leq t\leq G^* ,\quad \forall t\in G^0. \nonumber
\end{equation}
We further  define\quad $H^0=\{(\frac{\partial G}{\partial x})_i ,\quad i=1,2,.....,N \} $.
Let there exist some positive constant $W$ such that\quad $ {\left\vert t^0 \right \vert } \leq W ,\quad \forall \quad t^0\in H^0 $. So it is possible for very small $h_i,\forall i=1,2,...,N$,
\begin{equation}
  {\left\vert P_{lm} \right \vert }\leq 1+r_i, \quad \forall \quad i=l=m \quad l=1,2,.....,N, \qquad  \nonumber
\end{equation}
\begin{equation}
  {\left\vert P_{lm} \right \vert }\leq 1, \quad \forall\quad m=l+1,\quad i=l=1,2,.....,N-1,   \nonumber
\end{equation}
\begin{equation}
  {\left\vert P_{lm} \right \vert }\leq r_i, \quad \forall\quad i=l=m+1,\quad m=1,2,.....,N-2.   \nonumber
\end{equation}
Let \quad $ {\textbf{R}=[R_1,R_2,........,R_N]'_{1\times N}} $, denotes the row sum of the matrix $ \textbf{J}=(J_{lm})_{N\times N} $  where
\begin{equation}
  R_1=r_1+\frac{h^2_1}{24}(\alpha_1+\gamma_1)G_1+r_1\gamma_1\frac{h^3_1}{24}(\frac{\partial G}{\partial x})_i,  \quad l=i=1 , \qquad\qquad\qquad\qquad \qquad \nonumber
\end{equation}
\begin{equation}
  R_l=\frac{h^2_i}{24}(\alpha_i+\gamma_i+\beta_i)G_i+\frac{h^3_i}{24}(r_i\gamma_i-\beta_i)(\frac{\partial G}{\partial x})_i,\quad l=i=k, and \quad 2\leq k \leq N-1 , \qquad\quad \qquad\quad  \nonumber
\end{equation}
\begin{equation}
  R_N=1+\frac{h^2_N}{24}(\alpha_N+\beta_N)G_N-\frac{h^3_N}{24}\beta_N(\frac{\partial G}{\partial x})_N,\quad l=i=N . \qquad\qquad\qquad\qquad     \nonumber
\end{equation}

On neglecting the higher order terms i.e. $O(h^3_i)$ in $R_i$ then it is easy to see that $ \textbf{J} $ is irreducible \cite{Varga}. By the row sum criterion and for sufficiently small $h_i,\forall i=1,2,...,N$, $ \textbf{J} $ is monotone \cite{Henrici}. Thus $ \textbf{J}^{-1} $ exist and  $\textbf{J}^{-1}\geq 0 $. For the bound of $ \textbf{J} $, we define \cite{Volkov,Ahlberg}
\begin{equation}
  d_l(\textbf{J})=\left\vert J_{ll} \right \vert -\sum^{N}_{l\neq m}\left\vert J_{lm}\right\vert ,\quad l=1,2,....,N , \qquad\qquad    \nonumber
\end{equation}
where
\begin{equation}
   d_1(\textbf{J})=r_1+\frac{h^2_1}{24}(\alpha_1-\gamma_1)G_1-\frac{h^3_1}{6}r_1\gamma_1(\frac{\partial G}{\partial x})_1, \qquad\qquad\qquad\qquad\quad    \nonumber
\end{equation}
\begin{equation}
  d_l(\textbf{J})= \frac{h^2_i}{24}(\alpha_i-\beta_i-\gamma_i)G_i+\frac{h^3_i}{24}(\beta_i-r_i\gamma_i)(\frac{\partial G}{\partial x})_i,\quad l=i=k, and \quad 2\leq k\leq N-1 ,    \nonumber
\end{equation}
\begin{equation}
   d_N(\textbf{J})=1+\frac{h^2_N}{24}(\alpha_N-\beta_N)G_N+\frac{h^3_N}{24}\beta_N(\frac{\partial G}{\partial x})_N,\quad l=i=N .   \qquad  \nonumber
\end{equation}
We note  that higher order terms i.e. $O(h^3_i)$ in the above expressions are neglected. Let $d_l(\textbf{J})\geq0 ,\quad \forall\quad l$ \quad and
\begin{equation}
  d_*(\textbf{J})= \min_{\substack {1\leq l\leq N}}d_l(\textbf{J}).   \nonumber
\end{equation}
Then
\begin{equation}
 \|\textbf{J}^{-1}\|\leq \frac{1}{d_*(\textbf{J})}.
\end{equation}
Thus from (36) and (37), we have
\begin{equation}
  \|\textbf{E}\|\leq \frac{1}{d_*(\textbf{J})}\|\textbf{T}\|.
\end{equation}
It follows from (23) and (38) that  $\|\textbf{E}\|\rightarrow 0$ as $h_i\rightarrow 0$. Thus we conclude that method (11) converges and the order of the convergence of method (11) is at least cubic.

\section{Numerical Results}\label{sec4}
 To illustrate our method and demonstrate its computational efficiency, we considered some model problems. In each model problem, we took non-uniform step size $\emph{h}_i$. In Table 1 - Table 3, we have shown the maximum absolute error (MAY), computed for different values of \emph{N} and is defined as
 \begin{equation}
  MAY= \max_{\substack {1\leq{i}\leq {N}}}{\left\vert y(x_i)-y_i\right\vert }.       \nonumber
 \end{equation}

 The starting value of the step length $h_1$ is calculated by formula
 \begin{equation}
h_1=\begin{cases}
\frac{(b-a)(r-1)}{r^N-1} & \text{if  \quad $ r>1 $}\\
 \frac{(b-a)(1-r)}{1-r^N} & \text{if  \quad $ r<1 $ }         \nonumber
\end{cases}
 \end{equation}
 where $r=r_i ,\quad \forall\quad i=1,2,....,N $ in computation. The order of the convergence $(O_N)$ of the method (11) is estimated by the formula
 \begin{equation}
   (O_N)=\log_m(\frac{MAY_N}{MAY_{mN}}),   \nonumber
 \end{equation}
where m can be estimated by considering the ratio of $N's$.\\

We have used Newton-Raphson iteration method to solve the system of nonlinear equations  arisen from equation (11) or (12). All computations were performed on a Windows 7 Ultimate operating system in the GNU FORTRAN environment version 99 compiler  (2.95 of gcc) on  Intel Core i3-2330M 2.20 Ghz PC. The solutions are computed  on \emph{N} nodes and  iteration is continued until either the maximum difference between two successive iterates is less than $10^{-10}$ or the number of  iteration reached $10^5$. \\
\\
\textbf{Problem 1.} The first model problem is a linear problem \cite{Scott}  given by
\begin{multline}
  y''(x)=\frac{-3\epsilon}{(\epsilon+x)^2}y,\quad y(-0.1)=\frac{-0.1}{\sqrt{(\epsilon+0.01)}},\quad y(0.1)=\frac{0.1}{\sqrt{(\epsilon+0.01)}},\\ x\in[-0.1,0.1].      \nonumber
\end{multline}
The analytical solution  is $ y(x)=\frac{x}{\sqrt{(\epsilon+x^2)}}$. The MAY computed by method (11) for different values of \emph{N} and $\epsilon$ are presented in Table 1.\\
\textbf{Problem 2.} The second model problem is a nonlinear problem
\begin{equation}
  \epsilon y''(x)= \frac{3}{2}y^2 ,\quad y(0)=4,\quad y(1)= 1, \quad x\in[0,1].      \nonumber
\end{equation}
The analytical solution  is  $y(x)=\frac{4}{(1+\frac{x}{\sqrt(\varepsilon)})^2}$. The  MAY computed by method (11) for different values of \emph{N} are presented in Table 2.\\
\textbf{Problem 3.} The third model problem is a linear problem \cite{O'Riordan} given by
\begin{equation}
 \epsilon y''(x)=\frac{4}{(x+1)^4}(1+\sqrt{\epsilon}(x+1))y -f(x),\quad y(0)=2,\quad y(1)=-1,\quad x\in[0,1],      \nonumber
\end{equation}
where $f(x)$ is calculated so that $ y(x)=-\cos(\frac{4\pi x}{x+1})+\frac{3[\exp(\frac{-2\epsilon}{\sqrt{\epsilon}(x+1)})-\exp(\frac{-1}{\sqrt{\epsilon}})]}{1-\exp(\frac{-1}{\sqrt{\epsilon}})}$  is the analytical solution. The  MAY computed by method (11)for different values of \emph{N} and $\epsilon$ are presented in Table  3 .\\
\
\begin{table}[H]
\caption{Maximum absolute errors  (Problem 1).}
\vskip.4em
\setlength{\tabcolsep}{.10mm}
\renewcommand{\arraystretch}{1.25}
\begin{tabular}{|c|c|c|}
  \hline
  \multirow{2}{*}{} & {} & {Maximum absolute error }
       \\
       \setlength{\tabcolsep}{9.8 mm}
       ${\,r_1\,}$ & {\,\emph{N}\,}& \begin{tabular}{c c c c c}
             \hline
             % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
           ${\epsilon =1.0}$ \quad&\qquad${\epsilon =10^{-4}}$  \qquad \quad&\qquad ${\epsilon =10^{-6}}$ \qquad&\qquad ${\epsilon =10^{-8}}$\qquad &\quad\qquad ${\epsilon =10^{-10}}$ \\
             %\hline
           \end{tabular}
           \\  \cline {3-3}
           \hline
         \multirow{3}{*}{1.01}
         \setlength{\tabcolsep}{9.8 mm}
          & \,200 \, & \begin{tabular}{ c c c c c }
          .49479455(-2)\quad &\quad.18858910(-3)\quad&\quad.95963478(-5)\quad&\quad .26226044(-5)\quad&\quad .41723251(-6)\\
           \end{tabular}
           \\
           \setlength{\tabcolsep}{9.8 mm}
          & \,400 \, & \begin{tabular}{ c c c c c}
           .49697235(-2)\quad&\quad .15079975(-4)\quad&\quad .66757202(-5)\quad&\quad .71525574(-6)\quad&\quad .47683716(-6) \\
           \end{tabular}
           \\
           \setlength{\tabcolsep}{7.0mm}
          & \,800 \, & \begin{tabular}{ c c c c c }
             %\hline
             % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
            .50188862(-2)\quad&\quad.91791153(-5)\quad&\quad.67353249(-5)\quad&\quad.77486038(-6)\quad&\quad .10728836(-5) \\
           \end{tabular}
           \\
           \setlength{\tabcolsep}{7.0mm}
          & \,1600 \, & \begin{tabular}{ c c c c c }
             %\hline
             % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
            .50241202(-2)\quad&\quad.81062317(-5)\quad&\quad.31590462(-5)\quad&\quad.71525574(-6)\quad&\quad .23841858(-6) \\
           \end{tabular}
           \\ \cline {3-3}
           \hline
        \multirow{3}{*}{1.06}
         \setlength{\tabcolsep}{9.8 mm}
          & \,200 \, & \begin{tabular}{ c c c c c }
          .49364120(-2)\quad &\quad.29802322(-6)\quad&\quad.77486038(-6)\quad&\quad .11324883(-5)\quad&\quad .17881393(-6)\\
           \end{tabular}
           \\
           \setlength{\tabcolsep}{9.8 mm}
          & \,400 \, & \begin{tabular}{ c c c c c}
           .49356222(-2)\quad&\quad .35762787(-6)\quad&\quad .77486036(-6)\quad&\quad .11324883(-5)\quad&\quad .17881393(-6) \\
           \end{tabular}
           \\
           \setlength{\tabcolsep}{7.0mm}
          & \,800 \, & \begin{tabular}{ c c c c c }
             %\hline
             % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
            .49356222(-2)\quad&\quad .35762787(-6)\quad&\quad .77486036(-6)\quad&\quad .11324883(-5)\quad&\quad .17881393(-6) \\
           \end{tabular}
           \\ \cline {3-3}
           \hline
\end{tabular}
\end{table}
\
\begin{table}[H]
\caption{Maximum absolute errors with comparison to \cite{Pandey} (Problem 2).}
\vskip.4em
\setlength{\tabcolsep}{.10mm}
\renewcommand{\arraystretch}{1.25}
\begin{tabular}{|c|c|c|}
  \hline
  \multirow{2}{*}{} & {} & {Maximum absolute error }
       \\
       \setlength{\tabcolsep}{9.8 mm}
       ${\,r_1=1.08\,}$ & {\,\emph{N}\,}& \begin{tabular}{c c c c }
             \hline
             % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
         ${\epsilon =1.0}$  \qquad &\,\quad \qquad ${\epsilon =10.0}$ \qquad &\quad \qquad ${\epsilon =100.0}$\quad & \quad\qquad ${\epsilon =1000.0}$ \\
             %\hline
           \end{tabular}
           \\  \cline {3-3}
           \hline
           \multirow{3}{*}{}
         \setlength{\tabcolsep}{9.8 mm}
          & \,4 \, & \begin{tabular}{ c c c c  }
           .21700012(-2)\quad&\quad.81062317(-5)\quad&\quad .23841858(-6)\quad&\qquad.00000000 \\
           \end{tabular}
           \\
           \setlength{\tabcolsep}{9.8 mm}
          & \,8 \, & \begin{tabular}{ c c c c  }
           .75817108(-4)\quad&\quad.23841858(-6)\quad&\quad .23841858(-6)\quad&\quad.23841858(-6) \\
           \end{tabular}
           \\
           \setlength{\tabcolsep}{9.8 mm}
          & \,16 \, & \begin{tabular}{ c c c c  }
           .95367432(-6)\quad&\quad.23841858(-6)\quad&\quad.23841858(-6)\quad&\quad .23841858(-6) \\
           \end{tabular}
           \\
           \setlength{\tabcolsep}{9.8 mm}
         $method $ & \,32 \, & \begin{tabular}{ c c c c  }
           .23841858(-6)\quad&\quad.23841858(-6)\quad&\quad .23841858(-6)\quad&\quad .23841858(-6) \\
           \end{tabular}
           \\
           \setlength{\tabcolsep}{9.8 mm}
          (11) & \,64 \, & \begin{tabular}{ c c c c  }
           .23841858(-6)\quad&\quad.23841858(-6)\quad&\quad .23841858(-6)\quad&\quad .23841858(-6) \\
           \end{tabular}
           \\
           \setlength{\tabcolsep}{9.8 mm}
          & \,128 \, & \begin{tabular}{ c c c c  }
           .23841858(-6)\quad&\quad.23841858(-6)\quad&\quad .23841858(-6)\quad&\quad .23841858(-6) \\
           \end{tabular}
           \\
           \setlength{\tabcolsep}{9.8 mm}
          & \,256 \, & \begin{tabular}{ c c c c  }
          .11920929(-5)\quad& \quad .11920929(-5)\quad&\quad.11920929(-5)\quad&\quad .11920929(-5) \\
           \end{tabular}
           \\ \cline {3-3}
           \hline
           \multirow{3}{*}{}
         \setlength{\tabcolsep}{9.8 mm}
          & \,4 \, & \begin{tabular}{ c c c c  }
           .24600029(-1)\quad&\quad.12617111(-2)\quad&\quad .23126602(-4)\quad&\quad.23841858(-6) \\
           \end{tabular}
           \\
           \setlength{\tabcolsep}{9.8 mm}
          & \,8 \, & \begin{tabular}{ c c c c  }
           .60970783(-2)\quad&\quad.31232834(-3)\quad&\quad .52452087(-5)\quad&\quad.23841858(-6) \\
           \end{tabular}
           \\
           \setlength{\tabcolsep}{9.8 mm}
          & \,16 \, & \begin{tabular}{ c c c c  }
           .14939308(-2)\quad&\quad.87022781(-4)\quad&\quad.23841858(-6)\quad&\quad .23841858(-6) \\
           \end{tabular}
           \\
           \setlength{\tabcolsep}{9.8 mm}
      $ method \quad in\,$ & \,32 \, & \begin{tabular}{ c c c c  }
           .53775311(-3)\quad&\quad.31709671(-4)\quad&\quad .23841858(-6)\quad&\quad .23841858(-6) \\
           \end{tabular}
           \\
           \setlength{\tabcolsep}{9.8 mm}
         $\cite{Pandey}$ & \,64 \, & \begin{tabular}{ c c c c  }
           .36931038(-3)\quad&\quad.18835068(-4)\quad&\quad .23841858(-6)\quad&\quad .23841858(-6) \\
           \end{tabular}
           \\
           \setlength{\tabcolsep}{9.8 mm}
          & \,128 \, & \begin{tabular}{ c c c c  }
           .35333633(-3)\quad&\quad.19550323(-4)\quad&\quad .23841858(-6)\quad&\quad .23841858(-6) \\
           \end{tabular}
           \\
           \setlength{\tabcolsep}{9.8 mm}
          & \,256 \, & \begin{tabular}{ c c c c  }
          .35226345(-3)\quad& \quad.18119812(-4)\quad&\quad .23841858(-6)\quad&\quad .23841858(-6) \\
           \end{tabular}
           \\ \cline {3-3}
           \hline
\end{tabular}
\end{table}
\
\begin{table}[H]
\caption{Maximum absolute errors with comparison to \cite{Pandey} (Problem 3).}
\vskip.4em
\setlength{\tabcolsep}{.10mm}
\renewcommand{\arraystretch}{1.25}
\begin{tabular}{|c|c|c|}
  \hline
  \multirow{2}{*}{} & {} & {Maximum absolute error }
       \\
       \setlength{\tabcolsep}{9.8 mm}
       ${\,r_1=1.06\,}$ & ${\,\epsilon \,}$& \begin{tabular}{c c c c }
             \hline
             % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
         ${\emph{N}=100}$  \qquad &\,\quad \qquad ${ \emph{N}=500}$ \qquad & \qquad ${\emph{N}=1000}$\quad & \quad\qquad ${\emph{N}=1500}$ \\
             %\hline
           \end{tabular}
           \\  \cline {3-3}
           \hline
           \multirow{3}{*}{}
           \setlength{\tabcolsep}{9.8 mm}
          &$ \,1.0 \,$ & \begin{tabular}{ c c c c  }
          .32806396(-3)\quad&\quad.81896782(-4)\quad&\quad .8428968(-4)\quad&\quad .81658363(-4) \\
           \end{tabular}
           \\
         \setlength{\tabcolsep}{9.8 mm}
          & $ \,2^{-4} \,$ & \begin{tabular}{ c c c c  }
           .97751617(-5)\quad&\quad.11980534(-4)\quad&\quad .11861324(-4)\quad&\quad.12099743(-4) \\
           \end{tabular}
           \\
           \setlength{\tabcolsep}{9.8 mm}
          & $\,2^{-6} \,$ & \begin{tabular}{ c c c c  }
           .54240227(-5)\quad&\quad.56922436(-5)\quad&\quad .56624413(-5)\quad&\quad .54836273(-5) \\
           \end{tabular}
           \\
           \setlength{\tabcolsep}{9.8 mm}
          & $\,2^{-8} \,$ & \begin{tabular}{ c c c c  }
           .14781952(-4)\quad&\quad.14960766(-4)\quad&\quad .14930964(-4)\quad&\quad .14662743(-4) \\
           \end{tabular}
           \\
           \setlength{\tabcolsep}{9.8 mm}
          & $\,2^{-10} \,$ & \begin{tabular}{ c c c c  }
           .21040440(-4)\quad&\quad.20861626(-4)\quad&\quad .20891428(-4)\quad&\quad .20593405(-4) \\
           \end{tabular}
           \\
           \setlength{\tabcolsep}{9.8 mm}
         $method$ & $\,2^{-12} \,$ & \begin{tabular}{ c c c c  }
           .29563904(-4)\quad&\quad.29027462(-4)\quad&\quad .28967857(-4)\quad&\quad .28938055(-4) \\
           \end{tabular}
           \\
           \setlength{\tabcolsep}{9.8 mm}
         (11) &$ \,2^{-14} \,$ & \begin{tabular}{ c c c c  }
           .68480372(-1)\quad&\quad.65816164(-1)\quad&\quad .65812349(-1)\quad&\quad .65810800(-1) \\
           \end{tabular}
           \\
           \setlength{\tabcolsep}{9.8 mm}
          &$ \,2^{-16} \,$ & \begin{tabular}{ c c c c  }
          .97610575(0)\quad&\quad.94299066(0)\quad&\quad .94299060(0)\quad&\quad .94299066(0) \\
           \end{tabular}
           \\
           \setlength{\tabcolsep}{9.8 mm}
          & $\,2^{-18} \,$ & \begin{tabular}{ c c c c  }
          .27835369(-4)\quad&\quad.12516975(-5)\quad&\quad .26702881(-4)\quad&\quad .12516975(-5) \\
           \end{tabular}
           \\
           \setlength{\tabcolsep}{9.8 mm}
          & $\,2^{-20} \,$ & \begin{tabular}{ c c c c  }
          .22050738(-3)\quad&\quad.11920929(-5)\quad&\quad .11920929(-5)\quad&\quad .11920929(-5) \\
           \end{tabular}
           \\
           \setlength{\tabcolsep}{9.8 mm}
          & $\,2^{-22} \,$ & \begin{tabular}{ c c c c  }
          .26870966(-2)\quad&\quad.10728836(-5)\quad&\quad .10728836(-5)\quad&\quad .10728836(-5) \\
           \end{tabular}
           \\
           \setlength{\tabcolsep}{9.8 mm}
          &$ \,2^{-24} \,$ & \begin{tabular}{ c c c c  }
          .48569411(-1)\quad&\quad.11920929(-5)\quad&\quad .11920929(-5)\quad&\quad .11920929(-5) \\
           \end{tabular}
           \\ \cline {3-3}
           \hline
           \multirow{3}{*}{}
           \setlength{\tabcolsep}{9.8 mm}
          &$ \,1.0 \,$ & \begin{tabular}{ c c c c  }
          .38763934(-2)\quad&\quad.33957958(-2)\quad&\quad .33957660(-2)\quad&\quad .33954034(-2) \\
           \end{tabular}
           \\
         \setlength{\tabcolsep}{9.8 mm}
          & $ \,2^{-4} \,$ & \begin{tabular}{ c c c c  }
           .19446164(-2)\quad&\quad.19327551(-2)\quad&\quad .19329339(-2)\quad&\quad.19328594(-2) \\
           \end{tabular}
           \\
           \setlength{\tabcolsep}{9.8 mm}
          & $\,2^{-6} \,$ & \begin{tabular}{ c c c c  }
           .11941493(-2)\quad&\quad.11839569(-2)\quad&\quad .11841953(-2)\quad&\quad .11840165(-2) \\
           \end{tabular}
           \\
           \setlength{\tabcolsep}{9.8 mm}
          & $\,2^{-8} \,$ & \begin{tabular}{ c c c c  }
           .48494339(-3)\quad&\quad.48014522(-3)\quad&\quad .48032403(-3)\quad&\quad .48032403(-3) \\
           \end{tabular}
           \\
           \setlength{\tabcolsep}{9.8 mm}
          & $\,2^{-10} \,$ & \begin{tabular}{ c c c c  }
           .43356419(-3)\quad&\quad.36138296(-3)\quad&\quad .36138296(-3)\quad&\quad .36132336(-3) \\
           \end{tabular}
           \\
           \setlength{\tabcolsep}{9.8 mm}
         $ method\quad in\,$ & $\,2^{-12} \,$ & \begin{tabular}{ c c c c  }
           .49465895(-3)\quad&\quad.34594536(-3)\quad&\quad .34594536(-3)\quad&\quad .34594536(-3) \\
           \end{tabular}
           \\
           \setlength{\tabcolsep}{9.8 mm}
         $\cite{Pandey}$ &$ \,2^{-14} \,$ & \begin{tabular}{ c c c c  }
           .67013502(-3)\quad&\quad.33903122(-3)\quad&\quad .33909082(-3)\quad&\quad .33909082(-3) \\
           \end{tabular}
           \\
           \setlength{\tabcolsep}{9.8 mm}
          &$ \,2^{-16} \,$ & \begin{tabular}{ c c c c  }
          .11677146(-2)\quad&\quad.33509731(-3)\quad&\quad .33515692(-3)\quad&\quad .33521652(-3) \\
           \end{tabular}
           \\
           \setlength{\tabcolsep}{9.8 mm}
          & $\,2^{-18} \,$ & \begin{tabular}{ c c c c  }
          .27157217(-2)\quad&\quad.33360720(-3)\quad&\quad .33360720(-3)\quad&\quad .33372641(-3) \\
           \end{tabular}
           \\
           \setlength{\tabcolsep}{9.8 mm}
          & $\,2^{-20} \,$ & \begin{tabular}{ c c c c  }
          .79871528(-2)\quad&\quad.33271313(-3)\quad&\quad .33271313(-3)\quad&\quad .33271313(-3) \\
           \end{tabular}
           \\
           \setlength{\tabcolsep}{9.8 mm}
          & $\,2^{-22} \,$ & \begin{tabular}{ c c c c  }
          .25293380(-1)\quad&\quad.33223629(-3)\quad&\quad .33217669(-3)\quad&\quad .33223629(-3) \\
           \end{tabular}
           \\
           \setlength{\tabcolsep}{9.8 mm}
          &$ \,2^{-24} \,$ & \begin{tabular}{ c c c c  }
          .85839853(-1)\quad&\quad.33134222(-3)\quad&\quad .33128262(-3)\quad&\quad .33128262(-3) \\
           \end{tabular}
           \\ \cline {3-3}
           \hline
\end{tabular}
\end{table}
\
We have described a high order method for numerically solving two-point boundary value problems and several model problems considered to demonstrate the performance of the proposed method. Numerical result for examples 1 which is presented in table 1, for different values of $r_i$ show as $\epsilon $ decreases and $\emph{N} $ increases, for the non-uniform mesh size, maximum absolute error decreases but not substantial. The numerical results for examples 2, is maximum absolute error either decreases or remains same as there are change in \emph{N} but $\epsilon $ remains same. The results for examples 3, the maximum absolute error increases as $\epsilon $ decreases and \emph{N} increases. Over all method (11) is convergent and convergence of the method depends on choice of mesh ratio $ r_i$. The main advantage of the proposed high order method over method in \cite{Pandey} is its computational efficiency in considered model problems.
%It is an advantage of the exponential finite difference method over existing method \cite{Mohan}.


\section{ Conclusion}
A new non-uniform mesh size high order method to find the numerical solution of two point boundary value problems has been developed. At each nodal point $x = x_i,  i = 1, 2,....., N,$ we will obtain a system of algebraic equations given by (11). If the source function is f(x) then the system of equations from (11) is linear otherwise we will obtain nonlinear system of equations. It is obvious that special method required for some special problem where the solution is not regular and varies rapidly, so non-uniform mesh size method is natural choice. The new high order method produces good numerical approximate solutions for variety of model problems with non-uniform mesh i.e. $r_i\neq 1$. The numerical results of the model problems showed that the new high order method is computationally efficient. The rate of convergence of the present method is cubic.  The idea presented in this article leads to the possibility to develop non-uniform mesh size difference methods to solve third order and fourth order boundary value problems in ordinary differential equations. Works in these directions are in progress.

%\section{Bibliography}
%\ack The author is grateful to the anonymous reviewers and editor for their valuable suggestions, which substantially improved the standard of the paper.

\begin{thebibliography}{99}

\bibitem{Ascher} Ascher, U., Mattheij, R. M. M. AND Russell, R. D., Numerical Solution of Boundary Value Problems for Ordinary Differential Equations. Prentice-Hall, Englewood Cliffs, NJ (1988).

\bibitem{Collatz} Collatz, L ., Numerical Treatment of Differential Equations (3/ed.). Springer Verlag, Berlin,(1966).

\bibitem{Keller} Keller, H. B., Numerical Methods for Two Point Boundary Value Problems. Blaisdell Publ. Co.,New York, Waltham, Ma(1968).

\bibitem{Stoer} Stoer, J. and Bulirsch , R., Introduction to Numerical Analysis (2/e). Springer-Verlag, Berlin Heidelberg (1991).

\bibitem{Baxley} Baxley, J. V., Nonlinear Two Point Boundary Value Problems in Ordinary and Partial Differential Equations (Everrit, W.N. and Sleeman, B.D. Eds.), Springer-Verlag, New York (1981), 46-54.

\bibitem{Roos} Roos, H.G., Stynes, M. and Tobiska, L., Numerical Methods for Singularly Perturbed Differential Equations. Springer, New York (1996).

\bibitem{Hairer} Hairer, E., Nørsett, S. P.  and  Wanner, G., Solving Ordinary Differential Equations I Nonstiff Problems (Second Revised Edition). Springer-Verlag New York, Inc. New York, NY, USA (1993).

\bibitem{Jain} M.K. Jain, Numerical Solution of Difierential Equations, Wiley Eastern Limited, New Delhi, (1984).

\bibitem{Van}Van Niekerk, F.D., A nonlinear one step methods for initial value problems. Comp. Math. Appl. 13(1987), 367-371.

\bibitem{Pandey-1} Pandey, P. K., Solving Two Point Boundary Value Problems for Ordinary Differential Equations Using Exponential Finite Difference
      Method. Bol. Soc. Paran. Mat. (2014), **-**.

\bibitem{Pandey} Pandey, P. K. and Pandey B. D., Variable Mesh Size Exponential Finite Difference Method for the Numerical Solutions of Two Point Boundary Value Problems. Bol. Soc. Paran. Mat. (2014), **-**.

\bibitem{Varga} Varga, R. S., Matrix Iterative Analysis, Second Revised and Expanded Edition. Springer-Verlag, Heidelberg, (2000).

\bibitem{Golub}Golub, G. H. and Van Loan, C.F.,  Matrix Computations. John Hopkins University Press , Baltimore, USA (1989).

\bibitem{Henrici} Henrici, P., Discrete Variable Methods in Ordinary Differential Equations. John Wiley and Sons, New York, (1962).

\bibitem{Volkov} Volkov, Y. S. and Miroshnichenko, V. L. , Norm estimates for the inverses of matrices of monotone type and totally positive matrics. Siberian Mathematical Journal Vol. 50, No. 6,(2009), 982-987.

\bibitem{Ahlberg} Ahlberg, J. H. and Nilson, E. N. , Convergence properties of the spline fit. J. SIAM, 11, No. 1, (1963), 95-104.

\bibitem{Scott} Scott, M. R. and Watts, H. A., Computational solution of linear two point boundary value problems via orthonormalization. SIAM J. Numer. Anal. 14 (1977), 40-70.

\bibitem{O'Riordan} O'Riordan, E. and Stynes, M., A uniformely accurate finite element method for a singularly perturbed one-dimensional reaction-diffusion problem. Math. Comput.,47 (1986), 555-570.

\end{thebibliography}

\end{document}