\section{Introduction}
The field studying atomic gases at extremely low temperatures has seen riveting developments over the
past decade \cite{bloch2008}. This has pushed the field into a rather unique position where cold atoms
in various geometries can be used to do quantum simulation of widely used models from other fields including
condensed-matter physics \cite{lewenstein2007,esslinger2010}, particle and high-energy physics \cite{son2007,strings2008}, nuclear
physics \cite{zinner2013}, and even chemistry \cite{baranov2012}. A remarkably useful feature of cold atomic
gas experiments is the control over external trapping parameters that makes it possible to explore physics
in different dimensionalities \cite{bloch2008}. In particular, a number of impressive experiments have
studied various aspects of interacting bosons in one dimension \cite{moritz2003,stoferle2004,kinoshita2004,paredes2004,kinoshita2006,haller2009,haller2010}.
Using the atomic interaction resonances caused by one-dimensional
confinement \cite{olshanii1998} it has thus become possible to realize many interesting one-dimensional
systems including the famous hard-core bosonic Tonks-Girardeau gas \cite{tonks1936,girardeau1960} (see
Refs.~\cite{kinoshita2004,paredes2004}). Most recently, one-dimensional systems of interacting fermions have
been also realized in experiments \cite{pagano2014}.

In the last few years it has become possible to control the particle number in one-dimensional
experiments with great accuracy \cite{serwane2011}, allowing one to build a controllable few-body
system of fermions and study fermionization for strong interactions \cite{zurn2012},
pairing \cite{zurn2013}, impurity physics \cite{wenz2013}, a two-site Hubbard model \cite{murmann2015a},
and Heisenberg spin models \cite{murmann2015b}. These developments have sparked a great deal
of excitement in the community studying few-body physics and its relation to many-body phenomena.
Over the last decade a lot of new aspects of these one-dimensional systems have been covered by different authors.
This includes the physics of small trapped bosonic systems (single- and two-component)
\cite{zollner2005,deuret2007,tempfli2008,girardeau2011,brouzos2012,brouzos2014,wilson2014,zinner2014,garcia2015,dehk2015a,pietro2015},
details of few-fermion systems with two values of an internal degree of freedom (spin)
\cite{guan2009,yang2009,rubeni2012,brouzos2013,bugnion2013,gharashi2013,sowinski2013,gharashi2014,artem2014,cui2014,loft2015,lundmark2015,gharashi2015,unal2015,artem2015f} and various aspects of the transition from few- to many-body
physics \cite{girardeau2010,guan2010,astrak2013,volosniev2013,lindgren2014,deuret2014,volosniev2014,levinsen2014,sowinski2014,grining2015}. Furthermore, a number of studies have looked into mixed systems of bosons and
fermions and systems with particles of unequal mass
\cite{girardeau2004,girardeau2007,zollner2008,deuret2008,fang2011,harshman2012,garcia2013a,garcia2013b,harshman2014,campbell2014,garcia2014a,damico2014,mehta2014,garcia2014b,barf2015,artem2015b,dehk2015b,harshman2015,pecak2015,grass2015}.

In the present paper we are interested in studying one-dimensional
fermions or bosons with strong repulsive short-range interactions.
As has been recently discussed, in the strongly interacting limit the
system may be described by an effective Hamiltonian which has the
form of an (anisotropic) Heisenberg spin model \cite{volosniev2013,deuret2014}.
It was predicted theoretically that by tuning the interaction strength
from the weak (repulsive) to the strongly interacting regime for a
two-component Fermi system in one dimension one may arrive in the
ground state of an antiferromagnetic Heisenberg model \cite{volosniev2013}
and recently this was confirmed in experiments \cite{murmann2015b}.
However, the Heisenberg model obtained is non-trivial in the sense that
the exchange couplings are determined by the local geometry of the
trapping potential. Computing the local exchange coupling constants is
a formidable numerical task and thus several papers have discussed the
possibility to use various approximations to obtain these quantities.
A very neat approximation is the strong-coupling ansatz of Levinsen {\it et al.}
\cite{levinsen2014,pietro2015}
which allows one to get an extremely accurate set of exchange coefficients
for arbitrary system sizes in the case where the external confinement
is given by a harmonic oscillator. It has also been conjectured that for
smooth potentials, the well-known local density approximation (LDA) should give
a very accurate value for the local exchange coefficients \cite{deuret2014,levinsen2014}
and numerical results for up to six particles show that the deviations are
at most a few percent \cite{deuret2014}. This is a very reasonable
expectation and many studies in the past have used the LDA when
studying the properties of many-particle systems in one dimension
\cite{petrov2001,recati2003,astrak2004,tokatly2004,astrak2005,orso2007,liu2008,colo2008,rosi2015}.
However, this does not necessarily imply that the LDA works equally well
for smaller systems as large system sizes can sometimes average out some
of the finer details that the LDA may not capture.

In the present paper we test the performance of the LDA for strongly
interacting one-dimensional systems with up to six particles by
comparing it to exact calculations. We do this for several different
potential forms, including some experimentally relevant double-well
geometries. To the best
of our knowledge, no previous study has presented results for
the effective spin models with up to six particles in non-harmonic
confinement. We gauge the performance of the LDA against exact
results not only for statics (producing the local exchange constants
needed for the effective Hamiltonian) but also for dynamics. For
the latter case we study quantum state transfer using the spin models
that one obtains with the LDA and with exact calculations. Transfer of
quantum states in two-component spin systems is a delicate process that
depends sensitively on the local exchange couplings and thus provides a
difficult challenge for the LDA. The dynamical propagation
of information and correlations in one-dimensional setups with cold
atoms is a focal point of research at the
moment \cite{cheneau2012,fukuhara2013,hild2014,fukuhara2015}
and it is thus
theoretically important to have accurate models for these systems also
in the strongly interacting regime
where time-dependent exchange and non-equilibrium quantum magnetism can be
studied \cite{hild2014,trotzky2008,artem2015f}.

The paper is organized as follows. In Sec.~\ref{formal}
we outline the model and its assumptions, and then discuss how
to compute exchange coefficients exactly and within the LDA. Section~\ref{results} presents a
comparison of exact and LDA results for double-well potentials.
In Sec.~\ref{qst} we introduce the spin model picture of strongly
interacting 1D systems and we apply the local exchange coefficients
in the context of quantum state transfer to investigate some
consequences of the differences between different computational
schemes. Section~\ref{summary} contains a summary, discussions,
and outlook.

\section{Formalism}\label{formal}
We consider particles that are confined to move in one dimension (1D)
and assume that along the direction of motion there is a
trapping potential, $V(x)$, which is the same for all the
particles. The particles have short-range interactions
that we model by a Dirac delta-function and in turn the
Hamiltonian may be written in the following way
\begin{equation}
\label{hamiltonian}
H = \sum_{i=1}^N \left [ \frac{p_i^2}{2m} + V(x_i) \right ] + g_{1D} \sum_{i > j} \delta(x_i - x_j),
\end{equation}
where $p_i$ and $x_i$ are the momentum and coordinate operators of particle $i$, $m$ is the
mass of a particle and $V(x_i)$ is the trapping potential for the $i$th particle.
The interaction strength is parametrized by $g_{1D}$. In experimental setups the
1D confinement is achieved by applying a very tight transversal trap. It may
then be shown that $g_{1D}$ can be directly related to the scattering properties
of the non-trapped atoms and is a function of the low-energy three-dimensional
scattering length $a_{3D}$ and the transverse trapping length $a_\perp$ \cite{olshanii1998}.
What is very interesting is that one finds resonances where $g_{1D}$ diverges
due to the presence of the transverse confinement. This has been clearly demonstrated
in recent experiments \cite{zurn2012}. In particular, the regime where
$|g_{1D}|$ is very large is accessible experimentally \cite{murmann2015b}.
For simplicity we will use the notation $g=g_{1D}$ from now on since this
can give rise to no ambiguities in the present context.

In the present paper we will consider two-component Fermi systems with
$N = N_\uparrow + N_\downarrow$ atoms,
where $N_\uparrow$ and $N_\downarrow$ are number of atoms of components
with spin projection up and down, respectively. The zero-range two-body
interaction in the Hamiltonian will act only between pairs with opposite
spin projection. For pairs with the same spin projection the Pauli principle
requires antisymmetry upon exchange of the two particles. As the wave function
must also be continuous, we may infer that when two particles with identical
spin projections are close to each other the wave function has a continuous
first derivative. In turn, the Dirac delta-function two-body interaction
has no effect.
In experiment with atoms there are interactions between pairs of atoms with
identical spin projection, but they are highly suppressed due to their
short-range nature and thus can be safely neglected for our purposes.
The general Hamiltonian above takes this implicitly into account as
we ensure antisymmetry among like components in our $N$-body
wave functions. We note, however, that many of our results may
be easily transferred to two-component bosons with uniform
interactions, i.e. a single $g$ controlling the interactions
between both identical and different components \cite{volosniev2014}.

We now consider the strongly interacting regime where $1/g\to 0$
\cite{volosniev2013, volosniev2014, zinner2014}.
The most general eigenstate of the Hamiltonian in the
limit $1/g\to 0$ is \cite{volosniev2014}
\begin{equation}
\label{eigenstate}
\Psi = \sum_k a_k \theta(x_{P_k(1)}, \dotsc, x_{P_k(N)}) \Psi_0(x_1, \dotsc, x_N),
\end{equation}
where the sum is over the $N!$ permutations, $P_k$, of the coordinates, $a_k$ are the coefficients
that depend on the ordering of the particles, $\theta(x_1, \dotsc, x_N) = 1$, when $x_1 < x_2 < \dots < x_N$
and zero otherwise. $\Psi_0$ is a fully antisymmetric $N$-particle Slater determinant wave function constructed
from $N$ single-particle wave functions that are obtained by solving the
corresponding single-particle Schr{\"o}dinger equation with the potential $V(x)$. In
this paper we are interested in the lowest energy manifold of $N$-body states
which are obtained by taking the single-particle states with the $N$ lowest
energies.

The wave function $\Psi_0$ describes
the non-interacting $N$-fermion system with the energy $E_0$. It is
important to note that this $N$-body energy
is $M(N_\uparrow, N_\downarrow) = N! /(N_\uparrow!N_\downarrow!)$
times degenerate, and thus there is a manifold of $M(N_\uparrow,N_\downarrow)$
degenerate $N$-body states in the limit $1/g\to 0$ which is a quasi-degenerate
manifold at large but finite $g$ (which is where we will be working below as
we consider quantum state transfer). This degeneracy arises from the fact that
in the $1/g\to 0$ limit the particles become essentially impenetrable. Yet,
there are still $M(N_\uparrow,N_\downarrow)$ distinguishable ways that the particles may
be ordered on a line which all have the same energy. In practice one may
think of these various orderings of the spins along a line as a set of
basis states \cite{volosniev2013,volosniev2014,deuret2014}.

In the limit $1/g \to 0$ we can write the $N$-particle energy, $E$, to the linear
order in $1/g$ as~\cite{zinner2014, volosniev2013, volosniev2014}
\begin{equation}
E = E_0 - \frac{1}{g} \frac{\sum_{j=1}^{N-1} A_j \alpha_j}{\sum_{k=1}^{M(N_\uparrow, N_\downarrow)} a^2_k},
\end{equation}
where $A_j = \sum_{k>j} (a_{j} - a_{k})^2$. The important observation is that
there is an $a_k$-independent coefficient, $\alpha_j$, in this expression. It is
a geometric factor that depends solely on the total number of particles $N$ and
on the potential $V(x)$ and its single-particle
eigenstates. Remarkably, it does {\it not} depend on what system
one is considering as long as the interactions are strong, i.e. as long as we consider
the limit $1/g\to 0$. If one is able to compute $\alpha_j$ for a given $V(x)$, then
this may be used to study multi-component Fermi or Bose systems with all possible
combinations of internal components among the particles. The geometric factor
$\alpha_j$ can in a certain sense be thought of as the local exchange coupling
in the system. This interpretation is very useful when mapping the system onto
a spin model \cite{volosniev2013,volosniev2014,deuret2014,levinsen2014}. If
one thinks of the $N$ particles sitting on a line, then the index $j$ on
$\alpha_j$ corresponds to a pair of particles and the exchange coupling for
that pair is proportional to $\alpha_j$. Since strongly interacting 1D
systems are governed by exchange processes, the statics and dynamics of
such systems is essentially dictated by the $\alpha_j$ coefficients and
we would thus like to compute these in as general circumstances as
possible. Therefore, the coefficients $\alpha_j$ will therefore be the main focus
of our discussion.

An exact expression for $\alpha_j$ was first derived in \cite{volosniev2013}
\begin{equation}
\alpha_j =  \frac{\tfrac{\hbar^4}{m^2}\int \prod_{i=1}^{N}\mathrm{d}x_i \theta(x_1, \dotsc, x_N) \delta(x_1 - x_j)
\left(\partial \Psi_0\right)^{2}}{\langle\Psi\mid\Psi\rangle},
\label{eq:alphaExact}
\end{equation}
where the derivative in the integral denotes
\begin{equation}
\partial \Psi_0=\left[\frac{\partial \Psi_0}{\partial x_1}\right]_{x_1 = x_N},
\end{equation}
where one first takes the derivative of the $N$-body antisymmetric function
$\Psi_0$ with respect to $x_1$ and then subsequently sets $x_1=x_N$.
The exact expression is an $(N-1)$-dimensional integral and the numerical calculation of $\alpha_j$ is
by no means an easy task. It is therefore desirable to consider whether
appropriate approximations can be made to access these quantities also
for larger values of $N$. An important observation was made in Ref.~\cite{deuret2014}
where it was noticed that for $N\leq 6$, a local density approximation
can be used in computing $\alpha_j$ and this yields results that are
off by only a few percent for the benchmark case where $V(x)$ is a
simple harmonic oscillator potential. Using a highly accurate
ansatz wave function for the $N$-body problem, it later became
possible to get a highly accurate approximation to $\alpha_j$
for the harmonic trap for any $N$ given as a ratio of quadratic polynomials in $N$ \cite{levinsen2014}.

Here we are concerned with the question of how well the local density
approximation (LDA) does for different potentials. We must therefore
define and discuss how the LDA may be applied in the context
of strongly interacting particles in 1D. This discussion closely
parallels that of Ref.~\cite{deuret2014}. The main inspiration for the
LDA method in the present context comes from earlier work on the
Hubbard model using the Bethe ansatz where Ogata and Shiba have shown
that in the strongly interacting limit the spin and charge dynamics
decouple \cite{ogata1990}. The spin degrees of freedom may correspondingly
be described by a spin model of the Heisenberg type (we return to this
later on) with an exchange coupling that is proportional to the
third power of the density \cite{matveev2004,guan2007,matveev2008}.
This result is derived for a homogeneous system with periodic boundary
conditions which are the typical basic conditions needed to solve the
Bethe ansatz equations \cite{ogata1990}. In order to transfer these
results into the present context with non-homogeneous confinement
of the 1D system, Ref.~\cite{deuret2014} suggested to use the
density from the LDA in the expression for the exchange coupling
in the spin model. This yields the expression
\begin{equation}
\label{eq:alphaLDA}
\alpha^{(LDA)}_i = \frac{\hbar^4 \pi^2}{3m^2} n_{TF}^3(X_i).
\end{equation}
Here $n_{TF}(X_i)$ is the 1D Thomas-Fermi density~\cite{deuret2014,levinsen2014}
which is given by
\begin{equation}
n_{TF}(x) = \frac{1}{\pi\hbar} \sqrt{2m(\mu - V(x))},
\label{TFdensity}
\end{equation}
where $\mu$ is the chemical potential of the system. Ref.~\cite{deuret2014}
proposed to
calculate the Thomas-Fermi density
in the center-of-mass positions $X_i$ of the $i$th and $(i+1)$th particles.
The position $X_i$ can be
written as~\cite{deuret2014} $X_i = \frac{1}{2}\int \mathrm{d}x (\rho^{(i)}(x) + \rho^{(i+1)}(x))$,
where
\begin{equation}
\rho^{(i)}(x) = \frac{\int \mathrm{d}x_1 \dots \mathrm{d}x_N \theta(x_1, \dotsc, x_N) \delta(x - x_i)\Psi_0}{\langle\Psi\mid\Psi\rangle}.
\label{eq:density}
\end{equation}
Taking this position makes sense in light of the interpretation of
$\alpha$ as a local exchange coupling of a pair of particles and thus
by symmetry one should take the density at their common center-of-mass which
may be determined within the $N$-body system. One may argue that any
discrepancy between the exact value of $\alpha$ and the LDA approximation
can be alleviated by simply picking the right values of $X_i$. However,
finding the right $X_i$ appears as difficult as calculating the exact
$\alpha$ directly. Note also that for large $N$ finding $\rho^{(i)}$
is an equivalently difficult task and thus one needs some method of
determining densities accurately. Here we will ignore these problems
and be concerned with how well the LDA performs in cases where we can compute
both the exact and the LDA expression for $\alpha$ easily. We therefore
work with $N\leq 6$ here.
For a commonly used harmonic oscillator
trapping potential $V(x) = \frac{1}{2}m \omega^2 x^2$ the Thomas-Fermi
density can be written down in a simple form
\begin{equation}
\label{LDAHarmOsc}
n_{TF}(x) = \frac{1}{\pi \hbar}\sqrt{2m(N\hbar\omega - \frac{1}{2}m \omega^2 x^2)}.
\end{equation}
In this case the approximation is rather accurate
compared to the exact calculation and the relative error is not larger than $7$ percent and
decreases rapidly with increasing particle number $N$~\cite{deuret2014}.
However, for potentials of more complicated shape the calculation
of the center-of-mass positions becomes	much more difficult.

In the following we discuss the coefficients $\alpha_i$ and
$\alpha^{(LDA)}_i$, with $i = 1, \dotsc, N-1$ for different potential
profiles.
The numerical calculations of Eqs.~\eqref{eq:alphaExact} and~\eqref{eq:density} consist of two main parts.
First, we
construct the Slater determinant $\Psi_0$ with the eigenstates, $\psi_j(x_i)$, of a particle trapped in the potential $V(x_i)$.
We assume that our system is confined in a box potential with length $L$. To handle this numerically, we use the
normalized eigenstates
\begin{equation}
\phi_n(x_i) = \sqrt{\frac{2}{L}} \sin{\frac{\pi n x_i}{L}},
\end{equation}
where $n$ is an integer,
as an expansion basis for the functions
\begin{equation}
\psi_j(x_i) = \sum_n a^{(j)}_n \phi_n(x_i).
\end{equation}
To evaluate integrals in Eqs.~\eqref{eq:alphaExact}
and~\eqref{eq:density} we utilize a recursive many-dimensional integration
procedure with the one-dimensional integrals being calculated using a standard trapezoidal integration
routine. Even though more advanced integration methods of course are available they are not readily available
for the case of nested integrals (integration limits that depend on subsequent integrations).
We find that the trapezoidal integration
converges rapidly with the relative error for the values of the coefficients not exceeding $10^{-5}$. We find that
$50$ basis vectors $\phi_n$ and $50$ integration points are sufficient to achieve this accuracy in a reasonable time.
In the following sections we dismiss the superscript $(LDA)$ for
simplicity and make sure that it is clear from the context of the
discussion and in the figures whether we are discussing the exact
expressions for $\alpha$ or the LDA approximations.

\begin{figure}
\centering
\begin{minipage}{\linewidth}
\centering
\includegraphics[scale=0.65]{fig1a.eps}
\includegraphics[scale=0.65]{fig1b.eps}
\end{minipage}
\caption{The shapes of the potentials $V_1(x)$ from Eq.~\eqref{DW} (top panel) and $V_2(x)$ from Eq.~\eqref{AP} (bottom panel) for different values of parameters.}
\label{fig:potentials}
\end{figure}

\section{Comparison of exact and LDA coefficients}\label{results}
To perform a comparison
between the LDA results and the
exact solutions for the exchange coupling constants $\alpha_i$
we use two specific forms of a double well potential which have
flexibility to explore the similarities and differences between
the LDA and the exact calculation.
The first potential has the form
\begin{equation}
\label{DW}
V_{1}(x) = \frac{1}{2} k(|x| - b)^2,
\end{equation}
where $k$ and $b$ define the energy scale,
and the height of the central barrier. Notice that this
potential has a cusp (discontinuity of the first derivative)
at $x=0$. This is perfectly well allowed in the formalism as
this will not cause any troubling discontinuities in the
single-particle wave functions or their first derivatives.
The second potential is a
symmetric trap of the form~\cite{volosniev2014}
\begin{equation}
\label{AP}
V_{2}(x) = - V_{0} \sin^2 \left [ \frac{\pi}{2} \left ( \frac{2x}{L} + 1 \right )\right ] - u\sin^2 \left [ \pi \left (\frac{2x}{L} + 1\right)\right ],
\end{equation}
where the values of $V_0$ and $u$ control the shape of the potential.
Figure~\ref{fig:potentials} shows the shapes of the potentials we use in the
article for different values of the control parameters.

\begin{figure}
\begin{minipage}{\linewidth}
\includegraphics[scale=0.5]{fig2a.eps}
\includegraphics[scale=0.5]{fig2b.eps}
\includegraphics[scale=0.5]{fig2c.eps}
\end{minipage}\hfil
\caption{Geometric exchange coupling coefficients $\alpha_i$ for the
double well potential in Eq.~\eqref{DW} with $k=1$ as function of the
central barrier height $b$. The exact values based on the
expression in Eq.~\eqref{eq:alphaExact} are shown as solid lines, while the corresponding
approximate values based on the LDA formula in Eq.~\eqref{eq:alphaLDA} are the dashed lines.
The top panels shows $N=2,3$ and 4, the middle panel has $N=5$, and the bottom panel has
$N=6$. Note that the parity symmetry of Eq.~\eqref{DW} implies that $\alpha_{N-i} = \alpha_i$.
The dashed horizontal lines in the $N=4$ and $N=6$
cases correspond to the value $\alpha_1$ for $N=2$ and $N=3$,
respectively. This is the excepted asymptotic values
for large $b$ as discussed in the text.}
\label{fig:DWcomp}
\end{figure}

We note that we consider values in the interval $x \in [-L/2; L/2]$
where $L$ is the length of a box which confines the whole system.
This is required for our numerical procedure where we use
$L = 4 \pi$ in our calculations.
Energies are measured in units of $\epsilon=4\hbar^2/m L^2$,
with a factor of 4 coming from the box extending to $\pm L/2$
For the potential in Eq.~\eqref{DW} we set $k=1$ (in units of
$L$ and $\epsilon$) while for Eq.~\eqref{AP} we have used
$V_0=50$ in units of $\epsilon$. In deep enough potentials
the single-particle states do not feel the effect of this
'outside' boundary. However, for shallow potentials some of the states
might reside in the whole box and thus be influenced by the outer wall.
We will see how it affects the exchange constants in the subsections below.

\subsection{Results}
We now consider the coefficients $\alpha_i$ using Eqs.~\eqref{eq:alphaExact} and
~\eqref{eq:alphaLDA} for different values of the control parameters; the height of
the central barrier $b$ for the potential in Eq.~\eqref{DW} and the parameter $u$ for
the potential in Eq.~\eqref{AP}.
Figs.~\ref{fig:DWcomp} and~\ref{fig:APcomp} show both the LDA and exact
values of the coefficients $\alpha_i$ for different numbers of particles $N = 2-6$.
In this article we consider only spatially symmetric 1D potentials and so $\alpha_{N-i} = \alpha_i$.

In Fig.~\ref{fig:DWcomp} we show the results for the potential in Eq.~\eqref{DW}
comparing the exact results obtained from Eq.~\eqref{eq:alphaExact} to the LDA
results obtained with the formula in Eq.~\eqref{eq:alphaLDA} for particle numbers
$N=2-6$. Note that we do not need to plot all the $\alpha_i$ coefficients as the
parity invariance of the double well in Eq.~\eqref{DW} gives us the convenient
symmetry relation $\alpha_{N-i}=\alpha_i$. The results demonstrate that the LDA
does very well for small values of the central barrier height $b$ which should not
be too surprising as the potential for small $b$ resembles very much that of a
harmonic oscillator which was previously shown to be quite accurately described
within LDA \cite{deuret2014}.

One notices that for the even particle numbers the middle coefficient, i.e. $\alpha_1$,
$\alpha_2$, and $\alpha_3$ for $N=2$, $N=4$, and $N=6$, respectively, has similar behavior,
and that this behavior is not captured well by the LDA for $b\gtrsim 0.75$.
In fact, this middle coefficient, $\alpha_{N/2}$,
calculated via Eq.~\eqref{eq:alphaLDA} goes to zero much faster than the exact value.
This $\alpha_{N/2}$ pertains to the exchange coupling in the middle of the trap, i.e.
to the exchange taking place right at the central barrier in the potential. It is therefore
determined by the probability for particle tunneling through this central barrier, and
we clearly see the LDA fail to capture this effect for larger barriers, i.e. larger values
of $b$, where LDA underestimates exchange. This is connected to the fact that
LDA and in particular
the Thomas-Fermi density in Eq.~\eqref{TFdensity} is only well-defined in
between the classical turning points, where $\mu - V(x) > 0$. Hence for large barriers
where $\mu-V(x)\leq 0$, the Thomas-Fermi density goes to zero and thus in turn
the coefficient $\alpha_{N/2}$ goes to zero. For example, in the
case of the potential in Eq.~\eqref{DW} $\alpha^{(lda)}_{N/2} \equiv 0$ for $b > \sqrt{N}$ (when
$k=1$ as we have here). The termination values of $b=\sqrt{2}$, $b=2$, and $b=\sqrt{6}$ can be clearly
seen on the horizontal axis for $\alpha_1$ ($N=2$) and $\alpha_{2}$ ($N=4$) in the
top panel of Fig.~\ref{fig:DWcomp}, and for $\alpha_3$ ($N=6$) in the bottom panel of Fig.~\ref{fig:DWcomp}.

\begin{figure}
\begin{minipage}{\linewidth}
\includegraphics[scale=0.5]{fig3a.eps}
\includegraphics[scale=0.5]{fig3b.eps}
\includegraphics[scale=0.5]{fig3c.eps}
\end{minipage}\hfil
\caption{
Geometric exchange coupling coefficients $\alpha_i$ for the
smooth double well potential in Eq.~\eqref{AP} as function of the parameter $u$
that controls the central barrier. Here we have set $V_0=50$. The exact values based on the
expression in Eq.~\eqref{eq:alphaExact} are shown as solid lines, while the corresponding
approximate values based on the LDA formula in Eq.~\eqref{eq:alphaLDA} are shown as dashed lines.
The top panel shows $N=2,3$ and 4, the middle panel has $N=5$, and the bottom panel has
$N=6$. Note that the parity symmetry of Eq.~\eqref{AP} implies that $\alpha_{N-i} = \alpha_i$.}
\label{fig:APcomp}
\end{figure}

When the height of the barrier is large (large $b$) each of the two wells can be approximated by the harmonic oscillator
potential. We see that in this case the exact values of the coefficients becomes almost constant as function of the barrier height.
For even number of particles we intuitively expect that half of the particles will be located in each well. In the case with $N=4$
this would imply that $\alpha_1$ should go to a constant for large $b$, reaching the value corresponding to
the $b=0$ case with $N=2$ (two particles in a single harmonic well). This is clearly seen in the top panel
of Fig.~\ref{fig:DWcomp} where the horizontal dashed line marks the latter value. Likewise, for $N=6$, we expect a split into two three-body systems and thus that
both $\alpha_1$ and $\alpha_3$ should be approaching the $\alpha_1$ value for $b=0$ and $N=3$. Again this
is seen very nicely in the bottom panel of Fig.~\ref{fig:DWcomp} the first two coefficients for
$N=6$ approach the dashed horizontal line corresponding to $N=3$ and $b=0$ (notice that there
are different vertical scales on the three panels in Fig.~\ref{fig:DWcomp}).

For the spatially symmetric potentials considered here with an odd number of particles, we cannot use the same logic
of division of the particles into the two wells as the barrier grows large as for even particle numbers.
The overlap across the barrier of the single-particle wave functions will not vanish for odd particle numbers.
Hence, the coefficients $\alpha_i$ which depend solely on these single-particle wave functions will also not
vanish. We clearly see in the top and the middle panel of Fig.~\ref{fig:DWcomp} that for $b\gtrsim 2$
the LDA result does not do a good job in describing the exchange couplings. Again, this is caused by
the fast decrease of the Thomas-Fermi density on which the LDA result relies as the barrier increases.
This decrease of density clearly takes place much faster than seen in the exact results. We thus see
that odd-even effects can be considerable in the comparison of the LDA and the exact method.

In order to further explore the case with odd particle numbers we have tried to apply a small tilt
to the potential in order to explicitly break the parity invariance. This could for instance
be done by an additional term that is linear in $x$ in the potential~\eqref{DW}.
When this is done for an odd number of particles the values of the exchange couplings approach
the values of the couplings of a smaller system. For instance, for the system consisting of $N=5$
particles, two of the particles will be located in one of the wells, with the interacting coefficient approaching the value
of $\alpha_1$ for two particles in a harmonic trap, and the other three will occupy
the other well, with $\alpha_1 = \alpha_2$ approaching
those of a three-particle system in a harmonic trap.
At the same time we find the intuitive result, namely that
the interaction between these two subsystems will approach
zero as the barrier increases. We will not discuss the introduction
of slight symmetry-breaking terms any further here.
All in all, for the double well potential~\eqref{DW} the LDA approach provides a reliable approximation
to the exact values of the coefficients $\alpha_i$ for small central
barrier heights (small values of $b$) and then becomes poor for
larger barrier heights. An exception is found for even particles
numbers where the system splits into two equal size groups that
can be described as particles in two separate harmonic wells.
Here LDA does approach the exact results for the split system
asymptotically.

To further explore the robustness and/or failures of applying the LDA to
our problem system, we now change the potential into a different kind of
double well shape which has the form given in Eq.~\eqref{AP}. What one should
notice here is that the potential in Eq.~\eqref{AP} is in a
sense shallow, i.e. it does not increase to infinity around its edges.
Therefore the box potential that we have as a hard-wall boundary around the system
could be felt by the particles.
This is in fact the case for larger particle numbers as we will now discuss.

\begin{table}
\begin{center}
\begin{tabular}{| l | l | l |}
 \hline
 $N$&$\alpha_i^{(EXACT)} [\epsilon^2 L]$ & $\alpha_i^{(LDA)} [\epsilon^2 L]$\\ \hline \hline
 $2$ & $0.02487$ & $0.01328$ (47\%) \\ \hline
 $3$ & $0.06963$ & $0.04485$ (36\%) \\ \hline
 $4$ & $0.14921$ & $0.10630$ (29\%) \\ \hline
 $5$ & $0.27355$ & $0.20798$ (24\%) \\ \hline
 $6$ & $0.45260$ & $0.35911$ (21\%) \\
 \hline
\end{tabular}
\caption{Comparison of the exact calculation, $\alpha_{i}^{(EXACT)}$ (first column), to
that of the LDA, $\alpha_{i}^{(LDA)}$ (second column), for a potential consisting of a
flat box with infinite walls. Notice that in this case the exchange coefficients
are independent of $i$, i.e. $\alpha_i=\alpha$. The percentages in parenthesis in the
second column give the deviations of the LDA from the exact results. Notice that
the LDA always underestimates the exchange coupling coefficients. \label{tab:box}}
\end{center}
\end{table}

First we consider the cases with $N=2$ and $N=3$ as shown in the top panel
of Fig.~\ref{fig:APcomp}. Here we see that the qualitative behavior of the
$\alpha_i$ coefficients is similar to the potential of Eq.~\eqref{DW}
(shown in the top panel of Fig.~\ref{fig:DWcomp}) with the two-particle
case shown a steady decrease with the barrier parameter, $u$, while the
three-particle case first decreases and then has a slight increase again
at large values of $u$. This is analogous to the behavior in the top panel
of Fig.~\ref{fig:DWcomp}. We also notice that the LDA does a very good job
for the two- and three-particle systems for most values of $u$ in the
top panel of Fig.~\ref{fig:APcomp}. For the $N=3$ case this should be
compared and contrasted to the case in the top panel of Fig.~\ref{fig:DWcomp}
where the same $N=3$ case for the potential in Eq.~\eqref{DW} demonstrates
that the LDA deviates significantly from the exact result
for larger barrier heights (large $b$ values in that case).

For more particles ($N\geq 4$) the LDA results start to deviate significantly
from the exact results as we see in all three panels in Fig.~\ref{fig:APcomp}.
We see that already for four particles the coefficient $\alpha_i$
differs significantly in the LDA and exact calculations. This difference only increases
for larger particle numbers. The reason for this discrepancy can be traced back
to the form of the potential in Eq.~\eqref{AP} and its shallow nature as compared
to that in Eq.~\eqref{DW}. For the parameters we have chosen, the potential does
not become deep enough for five or six single-particle eigenstates to be confined
completely inside of the potential profile. This can be rephrased by stating that
not all the five or six lowest single-particle energy eigenvalues are in fact
negative. The most energetic single-particle wave functions that we require to
build the totally antisymmetric function, $\Psi_0$, are thus feeling the outside
of the potential. In particular, they are influenced by the box potential that
surrounds our system. While one may alleviate this problem by going to deeper
potentials (increasing $V_0$ in Eq.~\eqref{AP}) we choose $V_0=50$ in order
for us to study quantum state transfer in the last part of our paper and make
contact with recent results that utilize the
same parameters \cite{volosniev2014}.

To gain further insights into the influence of the box and the behavior of
LDA in this respect we may consider just the box potential on its own.
In Tab.~\ref{tab:box} we present the results of applying both the exact
formula (first column) and the LDA version (second column )to the box potential.
We clearly see a large discrepancy for the particle numbers we have studied.
The deviation of the LDA from the exact result is given in percentages in
the parenthesis in the second column of Tab.~\ref{tab:box}. By doing a
crude fit to the percentages we find that the deviations scale
approximately with the particle number as a power law $N^{-0.73}$. We thus
find a quite slow convergence and for smaller system sizes the deviations
can be significant.
We can see that as
the single-particle wave functions are distributed over the whole confining infinite square well potential
the local density approach becomes rather inadequate. The
large discrepancies in the values of the LDA and the exact results for
the $\alpha_i$ coefficients of the shallow potential in Eq.~\eqref{AP} for
particle numbers $N\geq 5$ can now
be better understood. They are in large part the result of the fact that the
outer box boundaries do become important for these particles numbers and
that the LDA does a poor job of describing particles in a box with a flat
bottom.
We do see
that the LDA approximation gives qualitatively similar results to the exact calculation even for $N > 4$, but
quantitatively the LDA may differ substantially. We have checked that as
one increases the depth of the potential in Eq.~\eqref{AP} one does indeed
see better agreement. However, the deviations we identified as common
for both potential profiles within the LDA remain.

\begin{figure}
\begin{minipage}{0.9\linewidth}
\includegraphics[scale=0.5]{fig4a.eps}
\includegraphics[scale=0.5]{fig4b.eps}
\includegraphics[scale=0.5]{fig4c.eps}
\includegraphics[scale=0.5]{fig4d.eps}
\end{minipage}
\caption{Quantum state transfer fidelity, $F$, for an $N=4$ two-component Fermi system with $J_i$ obtained from the
potential in Eq.~\eqref{AP} for $u=0$ (first row), $u=12.5$ (second row), $u=25$ (third row), and $u=50$ (fourth row)
where $u$ is in units of $\epsilon$. The other potential parameter is $V_0=50$ for all panels.
The (blue) solid line corresponds to the $J_i$ coefficients obtained from the exact calculation and
the (brown) dashed line corresponds to the LDA approach. Time is measured in units
of $\hbar/\epsilon$.}
\label{fig:APFidelityN=4}
\end{figure}

\section{Quantum transport properties}\label{qst}
In order to understand some potential effects that could be
implied by using the LDA instead of the exact solutions for the
exchange couplings, we now consider a dynamical protocol
that has recently been discussed in the context of strongly
interacting systems in 1D. As mentioned earlier, one
may indeed map these systems into effective spin models
\cite{volosniev2013,deuret2014}. This can be done for both
Fermi systems \cite{volosniev2013,deuret2014,levinsen2014}
and Bose systems \cite{volosniev2014,pietro2015}. To keep
the discussion concise, we will mainly consider the example of the
two-component Fermi system here. In that case, the spin
mapping is into the famous Heisenberg spin-$\tfrac{1}{2}$
model whose
Hamiltonian (up to a constant energy shift) is given by
\begin{equation}
 H_{s} =  \sum_{j=1}^{N-1}  J_j \bm S^j \cdot \bm S^{j+1},
 \label{eq:spinHam}
\end{equation}
where $\bm S^j = \frac{1}{2} \bm \sigma^j$ is a spin operator where $\bm \sigma^j = (\sigma_x^j, \sigma_y^j, \sigma_z^j)$
is the vector of the Pauli matrices. Here the nearest-neighbor
interaction coefficients are related very simply to the
local exchange coefficients and the coupling constant
as $J_{i} \equiv - \frac{\alpha_i}{g}$. This is another
justification for using the term 'local exchange
coupling' for the $\alpha_i$ above, i.e. that under
the spin mapping they appear as nearest-neighbor
couplings in what is equivalent to a spin chain
Hamiltonian. For Bose systems, one may have more
involved spin models as the coefficients for
$x$-, $y$-, and $z$- direction spin operators
are not necessarily the same (see Ref.~\cite{volosniev2014}
for further details). Below, we will make one
detour from the uniform Heisenberg spin model in
Eq.~\eqref{eq:spinHam} in order to consider the
so-called XX model which is a special case of
the Hamiltonian in Eq.~\eqref{eq:spinHam}
where all the terms with $z$-component operators
are eliminated.

\begin{figure}
\begin{minipage}{0.9\linewidth}
\includegraphics[scale=0.5]{fig5a.eps}
\includegraphics[scale=0.5]{fig5b.eps}
\includegraphics[scale=0.5]{fig5c.eps}
\includegraphics[scale=0.5]{fig5d.eps}
\end{minipage}\hfil
\caption{Same as Fig.~\ref{fig:APFidelityN=4} but for $N=5$.}
\label{fig:APFidelityN=5}
\end{figure}

Spin-$\tfrac{1}{2}$ chains have been proposed as
media for quantum state transfer about a decade ago
\cite{bose2003, bose2008, christandl2005}. The quantum
state transfer in these chains essentially corresponds
to flipping a single spin at one end of the chain
and then dynamically evolving the state such that one
may project onto the state at which the single spin
flip has reached the other end for all subsequent times.
The probability of finding the flipped spin at the other
end is known as the fidelity of the quantum state transfer
and if it reaches unity at some subsequent time we say
that the system supports perfect quantum state transfer
\cite{bose2008}. In the original proposal \cite{bose2003}
a spin chain with constant $J_i=J$ was considered and
it was shown that the fidelity could not reach unity
in this case. Later studies demonstrated that if the
coefficients are chosen as $J_i=J_0\sqrt{N(N-i)}$ with
$J_0$ some overall constant, then perfect state transfer
is possible (although only within the so-called XX model).
One would then have an ideal communications channel for quantum
information. However, it turns out to be exceedingly difficult
to produce a system that fulfills these requirements on
$J_i$. In Ref.~\cite{volosniev2014} Volosniev {\it et al.}
proposed strongly interacting 1D atomic systems as a possible
realization of perfect state transfer and found that the
potential in Eq.~\eqref{AP} can give rise to perfect state
transfer in the XX model for $u=12.5$. We will return to this below.
First we need to define the spin state space and the fidelity
that we consider.
For a spin chain with one impurity (a single spin that is flipped)
and $N-1$ majorities (in other words $N_\uparrow = 1$ and $N_\downarrow = N - 1$)
it is natural to use the basis of spin functions $\left \{ \mket{\uparrow\downarrow\downarrow\ldots\downarrow}, \mket{\downarrow\uparrow\downarrow\ldots\downarrow}
, \mket{\downarrow\ldots\downarrow\uparrow\downarrow}, \mket{\downarrow\ldots\downarrow\downarrow\uparrow}\right \}$. In this
basis we define the fidelity of the quantum state transfer which is a
function of time and is given by
\begin{equation}
 F(t) = \left | \mbra{\downarrow\ldots\downarrow\downarrow\uparrow} e^{-iH_{s}t/\hbar} \mket{\uparrow\downarrow\downarrow\ldots\downarrow} \right |^2.
 \label{eq:fidelity}
\end{equation}
The fidelity can be straightforwardly interpreted as the Hamiltonian
acting on the initial state on the right (with the single flipped
spin on the left edge of the systems) and then projection on the
final state on the left (with the single flipped spin on the
right edge). In practice, one expands the initial state on the
set of eigenstates of the Hamiltonian, then constructs the
state at all later times, and then projects onto the final
state which also has some expansion in terms of the
eigenstates of the Hamiltonian. Notice that the overall time
scale of the transfer depends on $\alpha_i$ and on $g$. In
the strict limit where $1/g=0$, there is no dynamics as
the particles are completely impenetrable and all orderings
are eigenstates (all of which are degenerate). However, in
the more realistic case where $g$ is large but finite, our
effective spin models work to linear order in $1/g$ and can
thus be used to study dynamics in the strongly interacting regime.
The timescale of transfer will then depend linearly on $g$.
One may think of the state transfer process as a set of
subsequent flips of pairs of spins along the chain. Each of
these local exchanges depend on $\alpha_i$, i.e. if $\alpha_i$
is large this happens fast and vice versa. We thus see that
large barriers will suppress the transfer as expected.
However, the process can depend rather delicately on the
actual values of the local exchanges as we will now demonstrate.
For instance, in the exact results we may see slow suppression
of the local exchange with barrier height as in Fig.~\ref{fig:DWcomp}
in regimes where LDA would give exactly zero (when the chemical
potential is below the barrier). This is one source of
error in the LDA that could carry into a dynamical protocol
like state transfer in a severe way.

\begin{figure}
\begin{minipage}{0.9\linewidth}
\includegraphics[scale=0.5]{fig6a.eps}
\includegraphics[scale=0.5]{fig6b.eps}
\includegraphics[scale=0.5]{fig6c.eps}
\includegraphics[scale=0.5]{fig6d.eps}
\end{minipage}\hfil
\caption{Same as Fig.~\ref{fig:APFidelityN=4} but for $N=6$.}
\label{fig:APFidelityN=6}
\end{figure}

Here we plot the fidelity of Eq.~\eqref{eq:fidelity} as a function of
time for different potentials of the form given in Eq.~\eqref{AP} for
the cases with $N=4$ (Fig.~\ref{fig:APFidelityN=4}), 5
(Fig.~\ref{fig:APFidelityN=5}), and 6 (Fig.~\ref{fig:APFidelityN=6}) particles,
respectively. All the three cases we present show the characteristic
oscillatory behavior of state transfer fidelities as function of
time \cite{volosniev2014}. However, one does indeed notice that
the oscillations become more prominent and more irregular as
$N$ increases (system size or chain length in the
spin chain language), and also as the control parameter for
the central barrier $u$ is increased.
Perfect state transfer
is never achieved for these parameters in a Heisenberg model of the type
in Eq.~\eqref{eq:spinHam}
without an external magnetic field. This is in accordance
with previous results~\cite{christandl2005,christandl2004,nikolopoulos2004}.
The calculations demonstrate that even when there are only small
differences in the exact and LDA results, the fidelity can show
large variations particularly for longer time intervals. This is
probably due to an accumulation of phase factors in the system over time that tends
to drive the exact and LDA results apart. However, we do see
some instances of very good agreement as for instance in the
second row of Fig.~\ref{fig:APFidelityN=4} where $u=12.5$
and $V_0=50$ in Eq.~\eqref{AP}.
Even in this case, noticeable
differences between exact and LDA results are seen,
but the overall agreement is very good.
For the larger
particle numbers in Figs.~\ref{fig:APFidelityN=5} and
\ref{fig:APFidelityN=6} we may even notice a clear
tendency for the LDA results to predict large fidelity
spikes that are either at different times as compared to
similar spikes using the exact result, or in some
instances the LDA shows spikes where the exact results
have none (see for instance the fourth row in
Fig.~\ref{fig:APFidelityN=5} or the third row in Fig.~\ref{fig:APFidelityN=6}).
We thus conclude that in most cases, the LDA results can
provide large deviations from the exact results in a
dynamical process such as quantum state transfer.

In closing this section, we want to consider how well
the LDA results do in the case where perfect state
transfer is achieved. This can be achieved with parameters
corresponding to the second row in Fig.~\ref{fig:APFidelityN=4}
in a setup where we can discard the $\bm S_z^j \bm S_z^{j+1}$
interaction in the Hamiltonian in Eq.~\eqref{eq:spinHam}
and thus reduce the problem to the Heisenberg XX model.
This can be realized by either using a tailored external magnetic
field applied to the system or in specific models with
strongly interacting two-component Bose systems \cite{volosniev2014}.
In Fig.~\ref{fig:APPerfTransferN=4} we compare the fidelities computed
using the
exact and the LDA values of the exchange coupling coefficients. We see that
the perfect transfer is indeed obtained for the values
calculated via Eq.~\eqref{eq:alphaExact}. However, the fidelity of
the spin transfer decreases with time for the coefficients
obtained via Eq.~\eqref{eq:alphaLDA}.
The relative differences in the $\alpha_i$ exchange coefficients
are no more than $2$ percent when comparing the exact to the LDA
method. However, we still see that
it affects the quantum transfer properties on
longer time scales significantly.

\begin{figure}
\centering
\includegraphics[scale=0.65]{fig7.eps}
\caption{Quantum state transfer fidelity, $F$, for an $N=4$ two-component Fermi
system with $J_i$ obtained from the
potential in Eq.~\eqref{AP} for $u=12.5$ and $V_0=50$ (in units of $\epsilon$)
for the case where the $z$-components in Eq.~\eqref{eq:spinHam}
can be discarded so that we are in the so-called Heisenberg XX spin model. In this case perfect state transfer
can be achieved \cite{volosniev2014}.
The (blue) solid line corresponds to the $J_i$ coefficients obtained from the exact calculation and
the (brown) dashed line corresponds to the LDA approach. Time is measured in units
of $\hbar/\epsilon$.}
\label{fig:APPerfTransferN=4}
\end{figure}

\section{Summary and discussion}\label{summary}
We have considered strongly interacting two-component systems in one dimension
held in place by an external confinement. Such systems have an effective Hamiltonian
that can be completely specified by computing a set of local exchange coefficients
which may be interpreted as nearest-neighbor spin exchange interactions when the
system is mapped onto a spin model Hamiltonian of the Heisenberg type. Computing
these exchange coefficients as accurately as possible is an important yet also
computationally difficult task. In the present paper we have explored two different
approaches. One is a 'brute force' calculation of a multi-dimensional integral which
gives the exact result (but is prohibitive for larger particle numbers) and the
other is an approach inspired by the local density approximation (LDA) that could in
principle be used to reduce the computational complexity. While previous studies
have shown that the local density approach can be accurate at the level of a
few percent for the case of a harmonic oscillator potential, we explore more
complicated geometries consisting of two instances of a double-well potential.

Our findings demonstrate that while the LDA does rather
well with potentials that resemble a harmonic oscillator, the
exchange couplings do not have the right qualitative and quantitative behavior
in the LDA when there are significant barriers as is typical of a double-well
potential. In particular, the LDA cannot capture the right quantum tunneling
processes across such barriers and may thus leave out important effects.
We have shown that for small systems this can lead to an underestimation of
exchange by the LDA for potentials which have a more interesting structure than the simple
single-well harmonic oscillator.
We expect this observation to have influence on experiments
with small system sizes in 1D optical lattices. There one typically also has an external
overall smooth potential (approximately harmonic in shape). Thus, the potential seen by
the particles is the superposition of harmonic trap and optical lattice, i.e. a 'smiling'
lattice potential.
This is a very structured potential and in the strongly interacting
regime one could be in dire straits with a simple LDA approach as compared to the exact exchange
couplings. In order to test out the LDA in a concrete physical process, we considered
quantum state transfer of single spin flips in systems with four, five, and six
particles. Here we found that while the LDA performs reasonably for some double-well
realizations there is amplification of the deviations of the LDA compared
to the exact results in the transfer fidelity that can be significant and lead to large
errors in both maximum fidelity values and the specific times at which these are
attained.

The current study has concentrated on small particle numbers where very accurate
results can be obtained for the local exchange coefficients using multi-dimensional
integration so that a comparison between exact and LDA results is possible for
arbitrary potentials. Incidentally, as we discussed in the introduction, the
system sizes used here are also of great current experimental interest in cold
atoms. One would, however, like to study how these results scale to larger particle
numbers. This most likely requires alternative approaches not only to the exact
formula but also to the LDA formula. In the current implementation it depends
on the density of the system which is not easy to compute for arbitrary potentials.
One may thus pursue an agenda of finding an alternative LDA method that obtains
the density by some other and computationally much simpler approach, and simultaneously
explore how to get a computational reduction of the exact formula.

The authors gratefully acknowledge discussions with and feedback from A.~G. Volosniev,
M. Valiente, D.~Petrosyan, N.~J.~S. Loft, A.~E. Thomsen, and L.~B. Kristensen. This
work was supported by the Danish Council for Independent Research DFF Natural Sciences
and the DFF Sapere Aude program.


\section{\textbf{Introduction}}
Nonlinear singular boundary value problems, commonly arises due to physical symmetry. Specifically, if a physical law is expressed by a Dirichlet problem $\bigtriangleup y(P)=f(P,y(P)),$ and one is interested in planar, cylindrical or spherical geometries, then it is led to the differential equation
\begin{eqnarray}
&&\label{1}-y''-\frac{\alpha}{x}y'=f(x,y),\quad 0 < x < 1,
\end{eqnarray}
where $\alpha= 0,1,2$, respectively \cite{Russell-Shampine-SIAM-JNA-1975}. In modern science various real life problems have been converted into a mathematical model similar to (\ref{1}), where $\alpha\geq0$, and several researchers have studied existence and uniqueness of the solution of (\ref{1}), e.g., (see \cite{Russell-Shampine-SIAM-JNA-1975, P.L.Chamber-JCP-1952, SC-DOVER-1967, J.B.Keller-JRMA-1956,Chawla-shivkumar-JCAM-1987, AKV-NA-2011}).

If the condition at the end point $(x=1)$ depends on some interior point of the interval $(0,1)$, then the following three point nonlinear SBVP arises
\begin{eqnarray}
\label{SBVP-1}&&-(x^{\alpha} y'(x))'= x^{\alpha}f(x,y) ,\quad 0<x<1,\\
\label{SBVP-2}&&y'(0)=0,\quad y(1)=\delta y(\eta),
\end{eqnarray}
where $f:I\times R\rightarrow R,~I=[0,1],~0<\eta<1,~\delta>0$ and $\alpha\geq1$. Different analytical techniques are available for solvability of three point nonlinear BVPs, for $\alpha=0, 1$ and $2$ (see \cite{FL-MJ-XL-CL-GL-NA-2008, J.R.L.Webb-NA-RWA-2012, Gregus-Neumann-Arscott-JLMS-1971, R.Ma-JMAA-2001, JN-NATMA-1997, JH-BK-CCT-PAMC-2005, AKV-MS-JOMC-2015}. The results of this paper may lead to some new developments in this area.

Two point BVPs are studied extensively when upper and lower solutions are well ordered. But when upper and lower solutions are in the reversed order (see \cite{AKV-MS-JOMC-2015, Zhang-SJMA-1995, Cabada-BVP-2011, Coster-Habets-Elsevier-2006, MC-CDC-PH-AMC-2001,CA-PA-SL-AMC-2001} and the reference there in) lot of explorations are still left unattended for both two point and multi point boundary value problems. In this paper we have tried to analyze few of these gaps related to monotone iterative techniques.

If $f(x,y)$ is continuous and Lipschitz continuous in its domain, we derive the following monotone iterative scheme from the nonlinear singular three point BVP (\ref{SBVP-1})--(\ref{SBVP-2}),
\begin{eqnarray}\label{MIS}
-y_{n+1}''(x)-\frac{\alpha}{x}y_{n+1}'(x)-\lambda y_{n+1}(x)=f(x,y_n)-\lambda y_n(x),\quad y_{n+1}'(0)=0,\quad y_{n+1}(1)=\delta y_{n+1}(\eta),
\end{eqnarray}
and prove that solution exists and belongs to the class $C[0,1]\cap~C^2[0,1]$.

We also define upper and lower solution which are our initial guesses for the above iterative scheme.
\begin{definition}
If the functions $u_0,v_0 \in C^2[0,1]$ are defined as
\begin{eqnarray}
\label{US-1}&& -(x^{\alpha}u_0'(x))'\geq x^{\alpha}f(x,u_0),\quad 0<x<1;\\
\label{US-2}&&~~u_0'(0)=0,\quad u_0(1)\geq \delta u_0(\eta),\\
\label{LS-1}&& -(x^{\alpha}v_0'(x))'\leq x^{\alpha}f(x,v_0),\quad 0<x<1;\\
\label{LS-2}&&~~v_0'(0)=0,\quad v_0(1)\leq \delta v_0(\eta),
\end{eqnarray}
 then $u_0$ and $v_0$ are called upper and lower solutions of the nonlinear three point SBVPs (\ref{SBVP-1})--(\ref{SBVP-2}), respectively.
\end{definition}
The purpose of this paper is to prove existence of unique solution for the class of nonlinear three point SBVPs. We observe that depending on the values of $\alpha$ we arrive at well ordered and reversed order cases. This classification, we deduce does not exist in the literature to the best of our knowledge.

This paper is organized in the following sections. Section \ref{TLC} we use Lommel's transformation to find out two linearly independent solutions in the terms of Bessel functions. Using these two linearly independent solutions Green's functions are constructed for different class of $\alpha$ (See Figure \ref{Figure}) in Section 3 and Section 4 states maximum  and anti-maximum  principles. Finally in Section 5 all these results are used to establish some new existence and uniqueness theorems. The sufficient conditions derived in this paper are verified for certain values which belongs to different classes of $\alpha$ in Section \ref{Example}.
\section{The Linear Case}\label{TLC}
The linear BVP corresponding to the nonlinear three point SBVPs (\ref{SBVP-1})--(\ref{SBVP-2}) is studied in this section. We consider the following inhomogeneous class of three point linear SBVPs,
\begin{eqnarray}
\label{Linear-SBVP}&&-(x^{\alpha} y'(x))'-\lambda x^{\alpha} y(x)=x^{\alpha} h(x),\quad 0<x<1,\\
&&\label{Linear-SBVP-B.C.-1}y'(0)=0,\quad y(1)=\delta y(\eta)+b,
\end{eqnarray}
where $h\in C(I)$ and $b$ is any constant. To solve the inhomogeneous system (\ref{Linear-SBVP})--(\ref{Linear-SBVP-B.C.-1}), we consider the corresponding homogeneous system
\begin{eqnarray}
\label{HP}&&-(x^{\alpha} y'(x))'-\lambda x^{\alpha} y(x)=0,\quad 0<x<1,\\
&&\label{HP-B.C.-1}y'(0)=0,\quad y(1)=\delta y(\eta).
\end{eqnarray}
Using Lommel's transformation (\S cf \cite{Chawla-shivkumar-JCAM-1987, A-Erdrlyi-McGraw-Hill-1953}) $z=\beta\zeta^{\gamma}$, $w= \zeta^{-a}v(\zeta)$, the standard Bessel's equation (\ref{2.4h}) is transformed into (\ref{2.5h})
\begin{eqnarray}
\label{2.4h}&&z^2\frac{d^2w}{dz^2}+z\frac{dw}{dz}+(z^2-\nu^2)w=0,\\
\label{2.5h}&&{\zeta}^2\frac{d^2v}{d{\zeta}^2}+{\zeta}(1-2a)\frac{dv}{d{\zeta}}+\left[(\beta \gamma \zeta^{\gamma})^2+(a^2-\nu^2\gamma^2)\right]v=0.
\end{eqnarray}
Now, by Lommel's Transformation the two linearly independent solutions of (\ref{2.5h}) are given by
\begin{eqnarray}
\label{2.6h}&&v_1(x)=\zeta^{a}w_1\left(\beta\zeta^{\gamma}\right),\quad v_2(x)= \zeta^{a}w_2\left(\beta\zeta^{\gamma}\right),
\end{eqnarray}
where $w_1(z)$ and $w_2(z)$ are two linearly independent solutions of Bessel's equation (\ref{2.4h}).\\
Now, if we set $\nu=a=\frac{1-\alpha}{2}$, $\gamma=1$, $\beta^2=\lambda$, then (\ref{2.5h}) reduces to (\ref{HP}) and hence we can obtained the two linearly independent solutions of (\ref{HP}) in terms of $w_1(z)$ and $w_2(z)$. A solution of (\ref{HP}) which is bounded in the neighborhood of the origin (except for a multiplicative constant) given by
$x^{\nu} J_{-\nu}\left(x \sqrt{\lambda}\right),~\mbox{if $\lambda>0$}~ \mathrm{and}~ x^{\nu} I_{-\nu}\left( x  \sqrt{|\lambda|}\right), ~\mbox{if $\lambda<0$}$.
\begin{note}
In this paper $J_{-\nu}$, $Y_{\nu}$ are Bessel functions of first and second kind and $I_{-\nu}$ and $K_{\nu}$ are Modified Bessel functions of first and second kind.
\end{note}
\section{Green's function} On the basis of sign of $\lambda$ and values of $\alpha$, we divide into the following cases.
\subsubsection{Case I: When $\lambda>0$ and $\alpha\notin \{1,3,5,\cdots\infty\}$.}
Suppose that
\begin{itemize}
\item [$(H_0):$]  $0<\lambda<j^{2}_{\nu,1}$,~~$0<\delta<1$,~~$\delta  \eta ^{\nu } J_{\nu }\left(\eta  \sqrt{\lambda }\right)-J_{\nu }\left(\sqrt{\lambda }\right)\geq 0$ and $J_{-\nu }\left(\sqrt{\lambda }\right)-\delta  \eta ^{\nu } J_{-\nu }\left(\eta  \sqrt{\lambda }\right)>0$,\\
      for $\alpha\in\bigcup_{n\in N} (4n-3,4n-1)$, and
\item [$(H_1):$]  $0<\lambda<j^{2}_{\nu,1}$,~~$\delta\geq1$,~~$\delta  \eta ^{\nu } J_{\nu }\left(\eta  \sqrt{\lambda }\right)-J_{\nu }\left(\sqrt{\lambda }\right)\leq 0$ and $J_{-\nu }\left(\sqrt{\lambda }\right)-\delta  \eta ^{\nu } J_{-\nu }\left(\eta  \sqrt{\lambda }\right)<0$,\\
     for $\alpha\in\bigcup_{n\in N} (4n-1,4n+1)$,
\end{itemize}
where $j^{2}_{\nu,1}$ is the first zero of $J_{\nu}(x)$. We can easily check the validation of $(H_0)$ and $(H_1)$ for $\alpha\in\bigcup_{n\in N} (4n-3,4n-1)$ and $\alpha\in\bigcup_{n\in N} (4n-1,4n+1)$ respectively.\\

Next two lemmas help us to define the sign of Green's function.

\pagebreak
\begin{figure}[H]
\includegraphics[width=16cm,height=21cm]{Flow-Chart-2}\\
\caption{Flow Chart}\label{Figure}
\end{figure}
\pagebreak

\begin{lemma}\label{LemmaJ-J-1} For $0<\lambda<\ j_{\nu,1}^2$, the Bessel functions of first kind $(J_{\nu}$ and $J_{-\nu})$ satisfy the following inequality
\begin{eqnarray*}
r^\nu \left(J_{-\nu}\left(s \sqrt{\lambda }\right)J_{\nu}\left(r \sqrt{\lambda }\right)-J_{\nu}\left(s \sqrt{\lambda }\right)J_{-\nu}\left(r \sqrt{\lambda }\right)\right)\geq 0,~~ 0< r\leq s\leq 1,
\end{eqnarray*}
where $\nu=\frac{1-\alpha}{2}$ and $\alpha\in\bigcup_{n\in N} (4n-3,4n-1)$.
\end{lemma}

\textbf{Proof.} Suppose
$$\widetilde{\Phi}(s,r)=r^\nu \left(J_{-\nu}\left(s \sqrt{\lambda }\right)J_{\nu}\left(r \sqrt{\lambda }\right)-J_{\nu}\left(s \sqrt{\lambda }\right)J_{-\nu}\left(r \sqrt{\lambda }\right)\right),$$ and let $s = s_0\in[0,1]$ be fixed. As
\begin{eqnarray*}
J_{-\nu}\left(s_0 \sqrt{\lambda }\right)J_{-1+\nu}\left(r \sqrt{\lambda }\right)+J_{\nu}\left(s_0 \sqrt{\lambda }\right)J_{1-\nu}\left(r \sqrt{\lambda }\right)&\leq&J_{-\nu}\left(s_0 \sqrt{\lambda }\right)J_{-1+\nu}\left(s_0 \sqrt{\lambda }\right)+J_{\nu}\left(s_0 \sqrt{\lambda }\right)J_{1-\nu}\left(x_0 \sqrt{\lambda }\right)\\
&=&\frac{2 \sin{\nu \pi}}{\pi s_0 \sqrt{\lambda}}\leq 0,
\end{eqnarray*}
for $r\leq s_0$. Now making use of the above inequalitiy, we can easily show that $\widetilde{\Phi}(s_0,t)$ is a non-increasing function of $r$. As $\widetilde{\Phi}(s_0,s_0)=0$ at $r=s_0$, which implies that $\widetilde{\Phi}(s_0,r)\geq 0,~~~\forall~~ r\leq s_0$. But as $s_0$ may have any value in $[0,1]$ therefore $\widetilde{\Phi}(s,r)\geq 0,~~\forall~~ 0< r\leq s\leq 1$.
\begin{lemma}\label{LemmaJ-J-2} For $0<\lambda<\ j_{\nu,1}^2$, the Bessel functions of first kind $(J_{\nu}$ and $J_{-\nu})$ satisfy the following inequality
\begin{eqnarray*}
r^\nu \left(J_{-\nu}\left(s \sqrt{\lambda }\right)J_{\nu}\left(r \sqrt{\lambda }\right)-J_{\nu}\left(s \sqrt{\lambda }\right)J_{-\nu}\left(r \sqrt{\lambda }\right)\right)\leq 0,~~ 0< r\leq s\leq 1,
\end{eqnarray*}
where $\nu=\frac{1-\alpha}{2}$ and $\alpha\in\bigcup_{n\in N} (4n-1,4n+1)$.
\end{lemma}

\textbf{Proof.} Proof follows the same analysis as we do in Lemma \ref{LemmaJ-J-1}, with the Bessel functions inequality
\begin{eqnarray*}
J_{-\nu}\left(s_0 \sqrt{\lambda }\right)J_{-1+\nu}\left(r \sqrt{\lambda }\right)+J_{\nu}\left(s_0 \sqrt{\lambda }\right)J_{1-\nu}\left(r \sqrt{\lambda }\right)&\geq&J_{-\nu}\left(s_0 \sqrt{\lambda }\right)J_{-1+\nu}\left(s_0 \sqrt{\lambda }\right)+J_{\nu}\left(s_0 \sqrt{\lambda }\right)J_{1-\nu}\left(s_0 \sqrt{\lambda }\right)\\
&=&\frac{2 \sin{\nu \pi}}{\pi s_0 \sqrt{\lambda}}\geq 0,
\end{eqnarray*}
for $r \leq s_0$.
\begin{lemma}\label{Lemma2.1}
For the linear three point SBVPs (\ref{HP})--(\ref{HP-B.C.-1}), where $\alpha\notin \{1,3,5,\cdots\infty\}$, the Green's function is given by
\begin{eqnarray*}
G(x,t)=\left\{
  \begin{array}{ll}
   \frac{\pi  \csc (\pi  \nu ) t^{\nu } x^{\nu } J_{-\nu }\left(x \sqrt{\lambda }\right) \left(\delta  \eta ^{\nu } \left(J_{\nu }\left(\eta  \sqrt{\lambda }\right) J_{-\nu }\left(t \sqrt{\lambda }\right)-J_{-\nu }\left(\eta  \sqrt{\lambda }\right) J_{\nu }\left(t \sqrt{\lambda }\right)\right)+\left(J_{-\nu }\left(\sqrt{\lambda }\right) J_{\nu }\left(t \sqrt{\lambda }\right)-J_{\nu }\left(\sqrt{\lambda }\right) J_{-\nu }\left(t \sqrt{\lambda }\right)\right)\right)}{2 \left(J_{-\nu }\left(\sqrt{\lambda }\right)- \delta  \eta ^{\nu } J_{-\nu }\left(\eta  \sqrt{\lambda }\right)\right)}, & 0 \leq x\leq t\leq \eta; \\
  \frac{1}{2} \pi  \csc (\pi  \nu ) t^{\nu } x^{\nu } J_{-\nu }\left(t \sqrt{\lambda }\right) \left(\frac{J_{-\nu }\left(x \sqrt{\lambda }\right) \left(\delta  \eta ^{\nu } J_{\nu }\left(\eta  \sqrt{\lambda }\right)-J_{\nu }\left(\sqrt{\lambda }\right)\right)}{J_{-\nu }\left(\sqrt{\lambda }\right)-\delta  \eta ^{\nu } J_{-\nu }\left(\eta  \sqrt{\lambda }\right)}+J_{\nu }\left(x \sqrt{\lambda }\right)\right), & t \leq x,~ t \leq \eta;\\
  \frac{\pi  \csc (\pi  \nu ) t^{\nu } x^{\nu } \left(J_{-\nu }\left(\sqrt{\lambda }\right) J_{\nu }\left(t \sqrt{\lambda }\right)-J_{\nu }\left(\sqrt{\lambda }\right) J_{-\nu }\left(t \sqrt{\lambda }\right)\right) J_{-\nu }\left(x \sqrt{\lambda }\right)}{2 \left(J_{-\nu }\left(\sqrt{\lambda }\right)- \delta  \eta ^{\nu } J_{-\nu }\left(\eta  \sqrt{\lambda }\right)\right)}, & x\leq t, ~\eta \leq t;\\
   \frac{1}{2} \pi  \csc (\pi  \nu ) t^{\nu } x^{\nu } \left(\frac{J_{-\nu }\left(x \sqrt{\lambda }\right) \left(\delta  \eta ^{\nu } J_{-\nu }\left(\eta  \sqrt{\lambda }\right) J_{\nu }\left(t \sqrt{\lambda }\right)-J_{\nu }\left(\sqrt{\lambda }\right) J_{-\nu }\left(t \sqrt{\lambda }\right)\right)}{J_{-\nu }\left(\sqrt{\lambda }\right)-\delta  \eta ^{\nu } J_{-\nu }\left(\eta  \sqrt{\lambda }\right)}+J_{-\nu }\left(t \sqrt{\lambda }\right) J_{\nu }\left(x \sqrt{\lambda }\right)\right), & \eta \leq t\leq x\leq 1.
  \end{array}
\right.
\end{eqnarray*}
If $H_0$ (or $H_1$) holds then $G(x,t)\leq0$ (or  $G(x,t)\geq0$).
\end{lemma}

\textbf{Proof.} We define Green's function as
\begin{eqnarray*}
G(x,t)=\left\{
  \begin{array}{ll}
   a_1~ x^{\nu } J_{-\nu }\left(x \sqrt{\lambda }\right)+a_2~ x^{\nu } J_{\nu }\left(x \sqrt{\lambda }\right) , & {0 \leq x\leq t\leq \eta;} \\
   a_3~ x^{\nu } J_{-\nu }\left(x \sqrt{\lambda }\right)+a_4~ x^{\nu } J_{\nu }\left(x \sqrt{\lambda }\right), & {t \leq x,~ t \leq \eta;} \\
   a_5~ x^{\nu } J_{-\nu }\left(x \sqrt{\lambda }\right)+a_6~ x^{\nu } J_{\nu }\left(x \sqrt{\lambda }\right), & {x \leq t,~ \eta \leq t;} \\
   a_7~ x^{\nu } J_{-\nu }\left(x \sqrt{\lambda }\right)+a_8~ x^{\nu } J_{\nu }\left(x \sqrt{\lambda }\right), & {\eta \leq t\leq x\leq 1}.
  \end{array}
\right.
\end{eqnarray*}
Using the following properties of the Green's function, for any $t \in [0,\eta]$ and boundary conditions, we have
\begin{eqnarray*}
&& a_1 t^{\nu}J_{-\nu}\left(t \sqrt{\lambda }\right)+ a_2 t^{\nu} J_{\nu}\left(t\sqrt{\lambda}\right)=a_3  t^{\nu}J_{-\nu}\left(t \sqrt{\lambda }\right)+ a_4 t^{\nu} J_{\nu}\left(t\sqrt{\lambda}\right),\\
&&-a_1 t^{\nu}\sqrt{\lambda }J_{1-\nu}\left(t \sqrt{\lambda }\right)+a_2 t^{\nu} \sqrt{\lambda} J_{\nu-1}\left(t\sqrt{\lambda}\right)+a_3 t^{\nu}\sqrt{\lambda }J_{1-\nu}\left(t \sqrt{\lambda }\right)-a_4 t^{\nu} \sqrt{\lambda} J_{\nu-1}\left(t\sqrt{\lambda}\right) = -\frac{1}{t^{1-2\nu}}
\end{eqnarray*}
and the following system of equations is derived
$$
\left(
  \begin{array}{cccc}
1 & 0 & -1 & 0\\
0 & 1 & 0 & -1\\
0 & 1 & 0 & 0\\
0 & 0 & J_{-\nu}\left( \sqrt{\lambda }\right)-\delta \eta^{\nu}J_{-\nu}\left(\eta \sqrt{\lambda }\right) & J_{\nu}\left(\sqrt{\lambda}\right)-\delta\eta^{\nu}J_{\nu}\left(\eta \sqrt{\lambda }\right)\\
  \end{array}
\right)
\left(
  \begin{array}{c}
    a_1 \\
    a_2 \\
    a_3 \\
    a_4 \\
  \end{array}
\right)=
\left(
  \begin{array}{c}
    \frac{\pi t^{\nu } J_{\nu }\left(t \sqrt{\lambda }\right)}{2 \sin{\nu \pi}} \\
    -\frac{\pi t^{\nu } J_{-\nu }\left(t \sqrt{\lambda }\right)}{2\sin{\nu\pi}} \\
      0 \\
      0 \\
 \end{array}
\right).$$
Solution of above system gives,
\begin{eqnarray*}
&&a_1=\frac{\pi  \csc (\pi  \nu ) t^{\nu } \left(J_{-\nu }\left(t \sqrt{\lambda }\right) \left(\delta  \eta ^{\nu } J_{\nu }\left(\eta  \sqrt{\lambda }\right)-J_{\nu }\left(\sqrt{\lambda }\right)\right)+J_{\nu }\left(t \sqrt{\lambda }\right) \left(J_{-\nu }\left(\sqrt{\lambda }\right)-\delta  \eta ^{\nu } J_{-\nu }\left(\eta  \sqrt{\lambda }\right)\right)\right)}{2 \left(J_{-\nu }\left(\sqrt{\lambda }\right)- \delta  \eta ^{\nu } J_{-\nu }\left(\eta  \sqrt{\lambda }\right)\right)},\\
&&a_2=0,\\
&&a_3=\frac{\pi  \csc (\pi  \nu ) t^{\nu } J_{-\nu }\left(t \sqrt{\lambda }\right) \left(\delta  \eta ^{\nu } J_{\nu }\left(\eta  \sqrt{\lambda }\right)-J_{\nu }\left(\sqrt{\lambda }\right)\right)}{2 \left(J_{-\nu }\left(\sqrt{\lambda }\right)- \delta  \eta ^{\nu } J_{-\nu }\left(\eta  \sqrt{\lambda }\right)\right)},\\
&&a_4=\frac{1}{2} \pi  \csc (\pi  \nu ) t^{\nu } J_{-\nu }\left(t \sqrt{\lambda }\right).
\end{eqnarray*}
Similarly for any $t \in [\eta,1]$, we have
\begin{eqnarray*}
&& a_5 t^{\nu}J_{-\nu}\left(t \sqrt{\lambda }\right)+ a_6 t^{\nu} J_{\nu}\left(t\sqrt{\lambda}\right)=a_7  t^{\nu}J_{-\nu}\left(t \sqrt{\lambda }\right)+ a_8 t^{\nu} J_{\nu}\left(t\sqrt{\lambda}\right),\\
&&-a_5 t^{\nu}\sqrt{\lambda }J_{1-\nu}\left(t \sqrt{\lambda }\right)+a_6 t^{\nu} \sqrt{\lambda} J_{\nu-1}\left(t\sqrt{\lambda}\right)+a_7 t^{\nu}\sqrt{\lambda }J_{1-\nu}\left(t \sqrt{\lambda }\right)-a_8 t^{\nu} \sqrt{\lambda} J_{\nu-1}\left(t\sqrt{\lambda}\right) = -\frac{1}{t^{1-2\nu}}.
\end{eqnarray*}
The above two equations and boundary conditions in $[\eta,1]$ gives
$$
\left(
  \begin{array}{cccc}
1 & 0 & -1 & 0\\
0 & 1 & 0 & -1\\
0 & 1 & 0 & 0\\
\delta \eta^{\nu}J_{-\nu}\left(\eta \sqrt{\lambda }\right) & \delta \eta^{\nu}J_{\nu}\left(\eta \sqrt{\lambda }\right) & -J_{-\nu}\left(\sqrt{\lambda }\right)  &  -J_{\nu}\left(\sqrt{\lambda }\right)\\
  \end{array}
\right)
\left(
  \begin{array}{c}
    a_5 \\
    a_6 \\
    a_7 \\
    a_8 \\
  \end{array}
\right)=
\left(
  \begin{array}{c}
    \frac{\pi t^{\nu } J_{\nu }\left(t \sqrt{\lambda }\right)}{2 \sin{\nu \pi}} \\
    -\frac{\pi t^{\nu } J_{-\nu }\left(t \sqrt{\lambda }\right)}{2\sin{\nu\pi}} \\
      0 \\
      0 \\
 \end{array}
\right).$$
By above four equation we have
\begin{eqnarray*}
&&a_5=\frac{\pi  \csc (\pi  \nu ) t^{\nu } \left(J_{-\nu }\left(\sqrt{\lambda }\right) J_{\nu }\left(t \sqrt{\lambda }\right)-J_{\nu }\left(\sqrt{\lambda }\right) J_{-\nu }\left(t \sqrt{\lambda }\right)\right)}{2 \left(J_{-\nu }\left(\sqrt{\lambda }\right)- \delta  \eta ^{\nu } J_{-\nu }\left(\eta  \sqrt{\lambda }\right)\right)},\\
&&a_6=0,\\
&&a_7=\frac{\pi  \csc (\pi  \nu ) t^{\nu } \left(\delta  \eta ^{\nu } J_{-\nu }\left(\eta  \sqrt{\lambda }\right) J_{\nu }\left(t \sqrt{\lambda }\right)-J_{\nu }\left(\sqrt{\lambda }\right) J_{-\nu }\left(t \sqrt{\lambda }\right)\right)}{2 \left(J_{-\nu }\left(\sqrt{\lambda }\right)- \delta  \eta ^{\nu } J_{-\nu }\left(\eta  \sqrt{\lambda }\right)\right)},\\
&&a_8=\frac{1}{2} \pi  \csc (\pi  \nu ) t^{\nu } J_{-\nu }\left(t \sqrt{\lambda }\right).
\end{eqnarray*}
This completes the construction of Green's function. Using $(H_0)$~(or $H_1$) and Lemma \ref{LemmaJ-J-1} (or Lemma \ref{LemmaJ-J-2}) we can easily verify that $G(x,t)\leq0$ (or  $G(x,t)\geq0$).
\subsubsection{Case II: When $\lambda>0$ and $\alpha\in \{1,3,5,\cdots\infty\}$.}
Suppose that
\begin{itemize}
\item [$(H_2):$]  $0<\lambda<y^{2}_{\nu,1}$,~~$0<\delta<1$,~~$\delta  \eta ^{\nu } Y_{\nu }\left(\eta  \sqrt{\lambda }\right)-Y_{\nu }\left(\sqrt{\lambda }\right)\geq 0$ and $J_{-\nu }\left(\sqrt{\lambda }\right)-\delta  \eta ^{\nu } J_{-\nu }\left(\eta  \sqrt{\lambda }\right)>0$, for~~$\alpha\in \{3,7,11,\cdots\}$, and
\item [$(H_3):$]  $0<\lambda<y^{2}_{\nu,1}$,~~$\delta\geq1$,~~$\delta  \eta ^{\nu } Y_{\nu }\left(\eta \sqrt{\lambda }\right)-Y_{\nu }\left(\sqrt{\lambda }\right)\leq 0$ and $J_{-\nu }\left(\sqrt{\lambda }\right)-\delta  \eta ^{\nu } J_{-\nu }\left(\eta \sqrt{\lambda }\right)<0$, for $\alpha\in\{1,5,9,\cdots\}$,

    where $y^{2}_{\nu,1}$ is the first zero of $Y_{\nu}(x)$.  It is easy to see that $(H_2)$ and $(H_3)$ can be satisfied for $\alpha\in \{3,7,11,\cdots\}$ and $\alpha\in \{1,5,9,\cdots\}$, respectively.
\end{itemize}

\begin{lemma}\label{LemmaJ-Y-1} For $0<\lambda<\ y_{\nu,1}^2$, the Bessel functions of first and second kind $(J_{-\nu}$ and $Y_{\nu})$ satisfy the following inequality
\begin{eqnarray*}
r^\nu\left(J_{-\nu}\left(s \sqrt{\lambda }\right)Y_{\nu}\left(r \sqrt{\lambda }\right)-Y_{\nu}\left(s \sqrt{\lambda }\right)J_{-\nu}\left(r \sqrt{\lambda }\right)\right)\geq 0,~~ 0 <r\leq s\leq 1,
\end{eqnarray*}
where $\nu=\frac{1-\alpha}{2}$ and $\alpha\in \{3,7,11,\cdots\}$.
\end{lemma}
\textbf{Proof.} Suppose
$$\widetilde{F}(s,r)=r^\nu\left(J_{-\nu}\left(s \sqrt{\lambda }\right)Y_{\nu}\left(r \sqrt{\lambda }\right)-Y_{\nu}\left(s \sqrt{\lambda }\right)J_{-\nu}\left(r \sqrt{\lambda }\right)\right),$$ and let $s = s_0\in[0,1]$ be fixed. Now as $J_{-\nu}\left(s_0 \sqrt{\lambda }\right)\geq J_{1-\nu}\left(r \sqrt{\lambda }\right)$ and $-Y_{\nu}\left(s_0 \sqrt{\lambda }\right) \geq Y_{-1+\nu}\left(r \sqrt{\lambda }\right)$ for $r \leq s_0$, when $0<\lambda\leq\ y_{v,1}^2$, where $\alpha\in \{3,7,11,\cdots\}$. So with the help of these inequalities we can easily show that $\widetilde{F}(s_0,r)$ is an non-increasing function of $r$. As $\widetilde{F}(s_0,s_0)=0$ at $r=s_0$, which implies that $\widetilde{F}(s_0,r)\geq 0,~~~\forall~~ r\leq s_0$. But as $s_0$ may take any value in $[0,1]$ therefore $\widetilde{F}(s,r)\geq 0,~~\forall~~ 0<r\leq s\leq 1$.

\begin{lemma}\label{LemmaJ-Y-2} For $0<\lambda<\ y_{\nu,1}^2$, the Bessel functions of first and second kind $(J_{-\nu}$ and $Y_{\nu})$ satisfy the following inequality
\begin{eqnarray*}
r^\nu\left(J_{-\nu}\left(s \sqrt{\lambda }\right)Y_{\nu}\left(r \sqrt{\lambda }\right)-Y_{\nu}\left(s \sqrt{\lambda }\right)J_{-\nu}\left(r \sqrt{\lambda }\right)\right)\leq 0,~~ 0< r\leq s\leq 1,
\end{eqnarray*}
where $\nu=\frac{1-\alpha}{2}$ and $\alpha\in\{1,5,9,\cdots\}$.
\end{lemma}

\textbf{Proof:} By using the inequalities $J_{-\nu}\left(s_0 \sqrt{\lambda }\right)\geq J_{1-\nu}\left(r \sqrt{\lambda }\right)$ and $Y_{\nu}\left(s_0 \sqrt{\lambda }\right) \geq -Y_{-1+\nu}\left(r \sqrt{\lambda }\right)$ for $r \leq s_0$, when $0<\lambda\leq\ y_{v,1}^2$, where $\alpha\in\{1,5,9,\cdots\}$, we can prove this lemma as we did in Lemma \ref{LemmaJ-Y-1}.

\begin{lemma}\label{Lemma2.2}
For the linear three point SBVPs (\ref{HP})--(\ref{HP-B.C.-1}), where $\alpha\in \{1,3,5,\cdots\infty\}$, the Green's function is given by
\begin{eqnarray*}
G(x,t)=\left\{
  \begin{array}{ll}
   \frac{\pi  \sec (\pi  \nu ) t^{\nu } x^{\nu } J_{-\nu }\left(x \sqrt{\lambda }\right) \left(J_{-\nu }\left(t \sqrt{\lambda }\right) \left(\delta  \eta ^{\nu } Y_{\nu }\left(\eta  \sqrt{\lambda }\right)-Y_{\nu }\left(\sqrt{\lambda }\right)\right)+Y_{\nu }\left(t \sqrt{\lambda }\right) \left(J_{-\nu }\left(\sqrt{\lambda }\right)-\delta  \eta ^{\nu } J_{-\nu }\left(\eta  \sqrt{\lambda }\right)\right)\right)}{2\left( J_{-\nu }\left(\sqrt{\lambda }\right)- \delta  \eta ^{\nu } J_{-\nu }\left(\eta  \sqrt{\lambda }\right)\right)}, & 0 \leq x\leq t\leq \eta; \\
 \frac{1}{2} \pi  \sec (\pi  \nu ) t^{\nu } x^{\nu } J_{-\nu }\left(t \sqrt{\lambda }\right) \left(\frac{J_{-\nu }\left(x \sqrt{\lambda }\right) \left(\delta  \eta ^{\nu } Y_{\nu }\left(\eta  \sqrt{\lambda }\right)-Y_{\nu }\left(\sqrt{\lambda }\right)\right)}{J_{-\nu }\left(\sqrt{\lambda }\right)-\delta  \eta ^{\nu } J_{-\nu }\left(\eta  \sqrt{\lambda }\right)}+Y_{\nu }\left(x \sqrt{\lambda }\right)\right), & t \leq x,~ t \leq \eta;\\
  \frac{\pi  \sec (\pi  \nu ) t^{\nu } x^{\nu } J_{-\nu }\left(x \sqrt{\lambda }\right) \left(J_{-\nu }\left(\sqrt{\lambda }\right) Y_{\nu }\left(t \sqrt{\lambda }\right)-Y_{\nu }\left(\sqrt{\lambda }\right) J_{-\nu }\left(t \sqrt{\lambda }\right)\right)}{2\left( J_{-\nu }\left(\sqrt{\lambda }\right)- \delta  \eta ^{\nu } J_{-\nu }\left(\eta  \sqrt{\lambda }\right)\right)}, & x\leq t, ~\eta \leq t;\\
  \frac{1}{2} \pi  \sec (\pi  \nu ) t^{\nu } x^{\nu } \left(\frac{J_{-\nu }\left(x \sqrt{\lambda }\right) \left(\delta  \eta ^{\nu } J_{-\nu }\left(\eta  \sqrt{\lambda }\right) Y_{\nu }\left(t \sqrt{\lambda }\right)-Y_{\nu }\left(\sqrt{\lambda }\right) J_{-\nu }\left(t \sqrt{\lambda }\right)\right)}{J_{-\nu }\left(\sqrt{\lambda }\right)-\delta  \eta ^{\nu } J_{-\nu }\left(\eta  \sqrt{\lambda }\right)}+J_{-\nu }\left(t \sqrt{\lambda }\right) Y_{\nu }\left(x \sqrt{\lambda }\right)\right), & \eta \leq t\leq x\leq 1.
  \end{array}
\right.
\end{eqnarray*}
If $H_2$ (or $H_3$) holds then $G(x,t)\leq0$ (or  $G(x,t)\geq0$).
\end{lemma}
\textbf{Proof:} The construction of Green's function follows the same analysis as we do in Lemma \ref{Lemma2.1}. Using the assumption $(H_2)$~(or $H_3$) and Lemma \ref{LemmaJ-Y-1} (or Lemma \ref{LemmaJ-Y-2}) we can easily verify that $G(x,t)\leq0$ (or  $G(x,t)\geq0$).

Now we state Lemmas \ref{Lemma2.3}, \ref{Lemma2.5} and \ref{Lemma2.6} and we omit proof.
\begin{lemma}\label{Lemma2.3}
If $y\in C^2(I)$ is a solution of inhomogeneous linear three point SBVPs (\ref{Linear-SBVP})--(\ref{Linear-SBVP-B.C.-1}), then
\begin{eqnarray}
\label{2.5}&&y(x)= \frac{b ~ x^{\nu}J_{-\nu}\left(x \sqrt{\lambda }\right)}{J_{-\nu }\left(\sqrt{\lambda }\right)-\delta  \eta ^{\nu } J_{-\nu }\left(\eta  \sqrt{\lambda }\right)}- {\int_0}^1{t^{\alpha}~ G(x,t) h(t)dt}.
\end{eqnarray}
\end{lemma}
\subsubsection{Case III: When $\lambda<0$.}
Suppose that
\begin{itemize}
\item [$(H'_0):$] $\delta>0$,~~$K_{\nu }\left(\sqrt{|\lambda| }\right)-\delta  \eta ^{\nu } K_{\nu }\left(\eta  \sqrt{|\lambda|  }\right)\leq 0$ and $I_{-\nu }\left(\sqrt{|\lambda|  }\right)-\delta  \eta ^{\nu } I_{-\nu }\left(\eta  \sqrt{|\lambda|  }\right)>0$ for~~$\alpha\in[1,\infty)$.
\end{itemize}

\begin{lemma}\label{LemmaI-K}
For $\lambda< 0$, the modified Bessel functions of first and second kind $(I_{-\nu}$ and $K_{\nu})$ satisfy the following inequality
\begin{eqnarray*}
r^\nu\left(K_{\nu }\left(s \sqrt{ |\lambda|  }\right) I_{-\nu }\left(r \sqrt{|\lambda|  }\right)-I_{-\nu }\left(s \sqrt{|\lambda|  }\right) K_{\nu }\left(r \sqrt{|\lambda|  }\right)\right)\leq 0,~~~ 0 <r\leq s\leq 1
\end{eqnarray*}
where $\nu=\frac{1-\alpha}{2}$ and $\alpha\in[1,\infty)$.
\end{lemma}
\textbf{Proof:}
Suppose
$$\widetilde{F_1}(s,r)=r^\nu\left(K_{\nu }\left(s \sqrt{ |\lambda|  }\right) I_{-\nu }\left(r \sqrt{|\lambda|  }\right)-I_{-\nu }\left(s \sqrt{|\lambda|  }\right) K_{\nu }\left(r \sqrt{|\lambda|  }\right)\right),$$ and further assume that $s = s_0\in[0,1]$ be fixed. Now we can easily show that the function $\widetilde{F_1}(s_0,r)$ will be non-decreasing for $r$ for all $\alpha\in[1,\infty)$. At $r=s_0$, $\widetilde{F_1}(s_0,s_0)=0$, i.e., $\widetilde{F_1}(s_0,r)\leq 0,~~\forall~~r\leq s_0$. But as $s_0$ may have any value in $[0,1]$ therefore $\widetilde{F}(s,r)\leq 0,~\forall~~ 0<r\leq s\leq 1$.
\begin{lemma}\label{Lemma2.5}
For the following linear three point SBVPs (\ref{HP})--(\ref{HP-B.C.-1}), where $\alpha\in [1,\infty)$, the Green's function is given by
\begin{eqnarray*}
G(x,t)=\left\{
  \begin{array}{ll}
   \frac{t^{\nu } x^{\nu } I_{-\nu }\left(x \sqrt{|\lambda|  }\right) \left(I_{-\nu }\left(t \sqrt{|\lambda|  }\right) \left(K_{\nu }\left(\sqrt{|\lambda|  }\right)-\delta  \eta ^{\nu } K_{\nu }\left(\eta  \sqrt{|\lambda|  }\right)\right)-K_{\nu }\left(t \sqrt{|\lambda|  }\right) \left(I_{-\nu }\left(\sqrt{|\lambda|  }\right)-\delta  \eta ^{\nu } I_{-\nu }\left(\eta  \sqrt{|\lambda|  }\right)\right)\right)}{I_{-\nu }\left(\sqrt{|\lambda|  }\right)-\delta  \eta ^{\nu } I_{-\nu }\left(\eta  \sqrt{|\lambda|  }\right)},~~~0 \leq x\leq t\leq \eta; \\
  \frac{t^{\nu } x^{\nu } I_{-\nu }\left(t \sqrt{|\lambda|  }\right) \left(I_{-\nu }\left(x \sqrt{|\lambda|  }\right) \left(K_{\nu }\left(\sqrt{|\lambda|  }\right)-\delta  \eta ^{\nu } K_{\nu }\left(\eta  \sqrt{|\lambda|  }\right)\right)-K_{\nu }\left(x \sqrt{|\lambda|  }\right) \left(I_{-\nu }\left(\sqrt{|\lambda|  }\right)-\delta  \eta ^{\nu } I_{-\nu }\left(\eta  \sqrt{|\lambda|  }\right)\right)\right)}{I_{-\nu }\left(\sqrt{|\lambda|  }\right)-\delta  \eta ^{\nu } I_{-\nu }\left(\eta  \sqrt{|\lambda|  }\right)},~~~ t \leq x,~ t \leq \eta;\\
  \frac{t^{\nu } x^{\nu } I_{-\nu }\left(x \sqrt{|\lambda|  }\right) \left(K_{\nu }\left(\sqrt{|\lambda|  }\right) I_{-\nu }\left(t \sqrt{|\lambda|  }\right)-I_{-\nu }\left(\sqrt{|\lambda|  }\right) K_{\nu }\left(t \sqrt{|\lambda|  }\right)\right)}{I_{-\nu }\left(\sqrt{|\lambda|  }\right)-\delta  \eta ^{\nu } I_{-\nu }\left(\eta  \sqrt{|\lambda|  }\right)}, ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ x\leq t, ~\eta \leq t;\\
    \frac{t^{\nu } x^{\nu } \left(I_{-\nu }\left(x \sqrt{|\lambda|  }\right) \left(K_{\nu }\left(\sqrt{|\lambda|  }\right) I_{-\nu }\left(t \sqrt{|\lambda|  }\right)-\delta  \eta ^{\nu } I_{-\nu }\left(\eta  \sqrt{|\lambda|  }\right) K_{\nu }\left(t \sqrt{|\lambda|  }\right)\right)-I_{-\nu }\left(t \sqrt{|\lambda|  }\right) K_{\nu }\left(x \sqrt{|\lambda|  }\right) \left(I_{-\nu }\left(\sqrt{|\lambda|  }\right)-\delta  \eta ^{\nu } I_{-\nu }\left(\eta  \sqrt{|\lambda|  }\right)\right)\right)}{I_{-\nu }\left(\sqrt{|\lambda|  }\right)-\delta  \eta ^{\nu } I_{-\nu }\left(\eta  \sqrt{|\lambda|  }\right)}, \\~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ \eta \leq t\leq x\leq 1,
  \end{array}
\right.
\end{eqnarray*}
and if $(H'_0)$ holds then $G(x,t)\leq0$.
\end{lemma}
\begin{lemma}\label{Lemma2.6}
If $y\in C^2(I)$ is a solution of inhomogeneous linear three point SBVPs (\ref{Linear-SBVP})-(\ref{Linear-SBVP-B.C.-1}) then
\begin{eqnarray*}
&&y(x)= \frac{b ~ x^{\nu}I_{-\nu}\left(x \sqrt{|\lambda|  }\right)}{I_{-\nu }\left(\sqrt{|\lambda|  }\right)-\delta  \eta ^{\nu } I_{-\nu }\left(\eta  \sqrt{|\lambda|  }\right)}- {\int_0}^1{t^{\alpha}~ G(x,t) h(t)dt}.
\end{eqnarray*}
\end{lemma}

\section{Maximum and anti-maximum principles for linear three point SBVPs}
The constant sign of Green's function results into the following anti-Maximum and maximum principles.
\begin{proposition}\label{Proposition2.1}
\textbf{$($Anti-maximum principle$)$}\\
Assume $\lambda>0$ and $(H_1)$ or $(H_3)$ holds, and $y\in C^2(I)$ satisfies
\begin{eqnarray*}
&&-(x^{\alpha} y'(x))'-\lambda x^{\alpha} y(x)\geq0,\quad 0<x<1,\\
&&y'(0)=0,\quad y(1)\geq \delta y(\eta).
\end{eqnarray*}
Then $y(x)\leq0$, $\forall x\in [0,1]$.
\end{proposition}
\begin{proposition}\label{Proposition2.2}
\textbf{$($Maximum principle$)$}
\begin{itemize}
\item[$(Max_1)$] Assume $\lambda>0$ and $(H_0)$ or $(H_2)$ holds, and $y\in C^2(I)$ satisfies
\begin{eqnarray*}
&&-(x^{\alpha} y'(x))'-\lambda x^{\alpha} y(x)\geq0,\quad 0<x<1,\\
&&y'(0)=0,\quad y(1)\geq \delta y(\eta).
\end{eqnarray*}
Then $y(x)\geq0$, $\forall x\in [0,1]$.
\item[$(Max_2)$]  Assume $\lambda<0$, $(H'_0)$ holds and $y\in C^2(I)$ satisfies
\begin{eqnarray*}
&&-(x^{\alpha} y'(x))'-\lambda x^{\alpha} y(x)\geq0,\quad 0<x<1,\\
&&y'(0)=0,\quad y(1)\geq \delta y(\eta).
\end{eqnarray*}
Then $y(x)\geq0$, $\forall x\in [0,1]$.
\end{itemize}
\end{proposition}
\section{\textbf{Existence results for nonlinear three point SBVP}}
The study of existence results for nonlinear three point SBVP is discussed in this section. On the basis of anti-maximum and maximum principles, we divide it into  the following two subsections.
\subsubsection{Reverse ordered upper $(u_0)$ and lower $(v_0)$ solutions $(u_0\leq v_0).$}
\begin{theorem}\label{Theorem3.1}
Assume that
\begin{itemize}
\item[$(R_1)$] there exist $u_0$, $v_0$ in $C^2[0,1]$ such that $u_0\leq v_0$, where $u_0$ satisfies (\ref{US-1})--(\ref{US-2}) and $v_0$ satisfies (\ref{LS-1})--(\ref{LS-2});
\item[$(R_2)$] the function $f: D \rightarrow R$ is continuous on $D :=\{(x,y) \in [0,1] \times R :u_0 \leq y \leq v_0\}$;
\item[$(R_3)$] there exist $M\geq 0$ such that for all $(x,y),(x,w)\in D$
\begin{eqnarray*}
y\leq w \Longrightarrow  f(x,w)- f(x,y) \leq M_1(w - y),
\end{eqnarray*}
\item[$(R_4)$] there exist a constant $\lambda>0$ such that $M_1-\lambda\leq0$ and $H_1$ or $H_3$ holds;
\end{itemize}
then the nonlinear three point SBVP (\ref{SBVP-1})--(\ref{SBVP-2}) has at least one solution in the region $D$. Sequence $\{u_n\}$ generated by equation (\ref{MIS}), with initial iterate $u_0$ converges monotonically (non-decreasing) and uniformly towards to the solution ${u(x)}$ of (\ref{SBVP-1})--(\ref{SBVP-2}). Similarly $v_0$ as an initial iterate leads to a non-increasing sequences $\{v_n\}$ converging to a solution $v(x)$. Any solution $z(x)$ in $D$ satisfies
\begin{eqnarray*}
 u(x)\leq z(x) \leq  v(x).
\end{eqnarray*}
\end{theorem}
\textbf{Proof:} It is easy to show that (see \cite{AKV-NA-2011})
\begin{eqnarray*}
u_{0}\leq u_{1}\leq u_{2}\leq  \ldots \leq u_{n}\leq u_{n+1}\leq  \ldots \leq v_{n+1}\leq v_{n}\leq \ldots \leq v_{2}\leq v_{1}\leq v_{0}.
\end{eqnarray*}
So the sequences $u_{n}$ is monotonically non-decreasing and bounded above by $v_0$, similarly $v_{n}$ is non-increasing, respectively and bounded below by $u_0$. Hence by Dini's theorem they converges uniformly. Let ${u}(x)=\displaystyle\lim_{n\to\infty}u_{n}(x) $ and ${v}(x)=\displaystyle\lim_{n\to\infty}v_{n}(x)$.

Now by using Lemma \ref{Lemma2.3}, the solution $u_{n+1}$ of (\ref{MIS}) is given by
\begin{eqnarray*}
u_{n+1}= \frac{b ~ x^{\nu}J_{-\nu}\left(x \sqrt{\lambda }\right)}{J_{-\nu }\left(\sqrt{\lambda }\right)-\delta  \eta ^{\nu } J_{-\nu }\left(\eta  \sqrt{\lambda }\right)}- {\int_0}^1{G(x,t)t^{\alpha} (f(t,u_{n})-\lambda u_{n})dt}.
\end{eqnarray*}
Now by Lebesgue's dominated convergence theorem, we get
\begin{eqnarray*}
&&{u}(x)= \frac{b ~ x^{\nu}J_{-\nu}\left(x \sqrt{\lambda }\right)}{J_{-\nu }\left(\sqrt{\lambda }\right)-\delta  \eta ^{\nu } J_{-\nu }\left(\eta  \sqrt{\lambda }\right)}- {\int_0}^1{G(x,t) t^{\alpha} (f(t,{u})-\lambda {u})dt}.
\end{eqnarray*}
Which is the solution of nonlinear SBVPs (\ref{SBVP-1})--(\ref{SBVP-2}). Similar equation we can define for the sequence of lower solution also. It is easy to see that ${u}(x)\leq z(x)\leq {v}(x)$.

\subsubsection{Well ordered upper $(u_0)$ and lower $(v_0)$ solutions $(v_0\leq u_0).$}
Based on the sign of $\lambda$, we prove two Existence theorems; Theorem \ref{Theorem3.2} and Theorem \ref{Theorem3.3}. The proof of these theorems are similar to the proof of Theorem \ref{Theorem3.1}.
\begin{theorem}\label{Theorem3.2}
Assume that
\begin{itemize}
\item[$(W_1)$] there exist  $u_0$, $v_0$ in $C^2[0,1]$ such that $v_0\leq u_0$, where $u_0$ satisfies (\ref{US-1})--(\ref{US-2}) and $v_0$ satisfies (\ref{LS-1})--(\ref{LS-2});
\item[$(W_2)$] the function $f: D_1 \rightarrow R$ is continuous on $D_1 :=\{(x,y) \in [0,1] \times R :v_0 \leq y \leq u_0\}$;
\item[$(W_3)$] there exist $M\geq 0$ such that for all $(x,y),(x,w)\in D_1$,
\begin{eqnarray*}
y\leq w \Longrightarrow  f(x,w)- f(x,y) \geq M_2(w - y),
\end{eqnarray*}
\item[$(W_4)$] there exist a constant $\lambda>0$ such that $M_2-\lambda\geq0$ and $H_0$ or $H_2$ holds;
\end{itemize}
then the nonlinear three point SBVP (\ref{SBVP-1})--(\ref{SBVP-2}) has at least one solution in the region $D_1$. Sequence $\{u_n\}$ generated by equation (\ref{MIS}), with initial iterate $u_0$ converges monotonically (non-increasing) and uniformly towards a solution ${u(x)}$ of (\ref{SBVP-1})--(\ref{SBVP-2}). Similarly $v_0$ as an initial iterate leads to a non-decreasing sequences $\{v_n\}$ converging to a solution $v(x)$. Any solution $z(x)$ in $D_1$ satisfies
\begin{eqnarray*}
v(x)\leq  z(x) \leq u(x).
\end{eqnarray*}
\end{theorem}
\begin{theorem}\label{Theorem3.3}
Assume that
\begin{itemize}
\item[$(W'_1)$] there exist  $u_0$, $v_0$  in $C^2[0,1]$ such that $v_0\leq u_0$, where $u_0$ satisfies (\ref{US-1})--(\ref{US-2}) and $v_0$ satisfies (\ref{LS-1})--(\ref{LS-2});
\item[$(W'_2)$] the function $f: D_2 \rightarrow R$ is continuous on $D_2:=\{(x,y) \in [0,1] \times R :v_0 \leq y \leq u_0\}$;
\item[$(W'_3)$] there exist $ M_3 \geq 0$ such that for all $(x,\widetilde{y}),(x,\widetilde{w})\in D_2$
\begin{eqnarray*}
\widetilde{y}\leq \widetilde{w} \Longrightarrow  f(x,\widetilde{w})- f(x,\widetilde{y}) \geq  -M(\widetilde{w} - \widetilde{y})
\end{eqnarray*}
\item[$(W'_4)$] there exist a constant $\lambda<0$ such that $ M_3+\lambda \leq 0$ and $(H_0')$ holds;
\end{itemize}
then the nonlinear three point SBVP (\ref{SBVP-1})--(\ref{SBVP-2}) has at least one solution in the region $D_2$. Sequences $\{u_n\}$ generated by equation (\ref{MIS}), with initial iterate $u_0$ converges monotonically (non-increasing) and uniformly towards a solution ${u(x)}$ of (\ref{SBVP-1})--(\ref{SBVP-2}). Similarly $v_0$ as an initial iterate leads to a non-decreasing sequences $\{v_n\}$ converging to a solution $v(x)$. Any solution $ z(x)$ in $D_2$ satisfies
\begin{eqnarray*}
v(x)\leq z(x) \leq u(x).
\end{eqnarray*}
\end{theorem}

\subsection{Uniqueness of nonlinear three point SBVP}
\begin{theorem}
Let $f(x,y)$ be continuous on $D$ (or $D_1$ or $D_2$) and there exist a constant $M_{\lambda}\geq0$ such that
\begin{eqnarray*}
&& f(x,u)- f(x,v) \leq M_{\lambda}(u - v),
\end{eqnarray*}
and $M_{\lambda}<\lambda_1$, where $\lambda_1 \leq \min{\left\{j^{2}_{\nu,1},y^{2}_{\nu,1},i^{2}_{-\nu,1},k^{2}_{\nu,1}\right\}}.$ Then the nonlinear three point SBVP (\ref{SBVP-1})--(\ref{SBVP-2}) has a unique solution.
\end{theorem}
\textbf{Proof:} Suppose $u(x)$ and $v(x)$ be any two solutions of (\ref{SBVP-1})--(\ref{SBVP-2}) then we get
\begin{eqnarray*}
&&-(x^{\alpha} (u-v)')'= x^{\alpha}[f(x,u)-f(x,v)],\\
&&(u-v)'(0)=0,~~~~(u-v)(1)= \delta (u-v)(\eta),
\end{eqnarray*}
which gives
\begin{eqnarray*}
&&-(x^{\alpha} (u-v)')'- M_{\lambda} x^{\alpha}(u-v)\leq 0,\\
&& (u-v)'(0)=0,~~~~(u-v)(1)= \delta (u-v)(\eta).
\end{eqnarray*}
Since by the Maximum and Anti-maximum principles $(b=0)$, whenever $M_{\lambda}<\lambda_1$, we get $u-v\leq 0$ or $u-v\geq 0$ (i.e., $u\leq v$ or $u\geq v$) for different class of $\alpha$. Similarly by changing the role of $u$ and $v$, we get $u\geq v$ or $u\leq v$. Hence $u \equiv v$. Therefor the solution of the (\ref{SBVP-1})--(\ref{SBVP-2}) is unique.

\section{Numerical illustrations}\label{Example}
We present here some numerical examples to validate our existence results which is derived in the Theorem \ref{Theorem3.1}, Theorem \ref{Theorem3.2} and Theorem \ref{Theorem3.3}.
\subsection{Reverse ordered upper and lower solution}
The following examples validate the result of Theorem \ref{Theorem3.1}, and gives a range of $\lambda$ for that we can generate two monotone sequences which converge to the solution of nonlinear SBVP.
\begin{example}
Consider the nonlinear three point SBVP
\begin{eqnarray}
\label{Ex4.1.1}&& -y''(x)-\frac{\alpha}{x} y'(x)=\frac{\alpha~(e^{y}-1)-x}{4},\\
\label{Ex4.2.1}&&~~y'(0)=0,~~~~~~ y(1)=3 y\left(\frac{1}{7}\right).
\end{eqnarray}
\end{example}
Where $f(x,y)=\frac{\alpha~(e^{y}-1)-x}{4},~~\delta=3,~\eta =\frac{1}{7}$ and $\alpha$ satisfies $(H_1)$ or $(H_3)$. In this problem we choose lower and upper solutions as $v_0= 1$ and $u_0=-1$, where $v_0\geq u_0$, i.e., it is reverse ordered case. Here nonlinear term satisfies all assumptions for Theorem \ref{Theorem3.1} and Lipschitz constant is $M_1=\frac{e\alpha}{4}$. Now we can find out a range for $\lambda\in (\frac{e\alpha}{4},~~j_{\nu,1}^2)$ such that $(H_1)$ or $(H_3)$ are true.
\subsection{Well ordered upper and lower solution}
The following examples validate the results of Theorem \ref{Theorem3.2} and Theorem \ref{Theorem3.3}. On the basis of sign of $``\lambda"$, we divide this subsection into the following two parts.
\subsubsection{When $\lambda >0$}
\begin{example}
Consider the nonlinear three point SBVP
\begin{eqnarray}
\label{Ex4.1.2}&& -y''(x)-\frac{\alpha}{x} y'(x)=de^{y},\\
\label{Ex4.2.2}&&~~y'(0)=0,~~~~~~ y(1)=\frac{1}{2} y\left(\frac{1}{3}\right).
\end{eqnarray}
\end{example}
where $f(x,y)=de^{y},~~d=\frac{2(1+\alpha)e^{-\frac{2}{3}}}{3}~~\delta=\frac{1}{2},~\eta =\frac{1}{3}$ and $\alpha$ satisfies $(H_0)$ or $(H_2)$. In this problem we choose lower and upper solutions as $v_0= 0$ and $u_0=\frac{2-x^2}{3}$, where $v_0\leq u_0$, i.e., it is well ordered case. The nonlinear term satisfies all assumptions for Theorem \ref{Theorem3.2} and Lipschitz constant is $M_2=d$. Now we can find out a range for $\lambda>0$ such that the conditions $M_2-\lambda \geq 0$,~$(H_0)$ or $(H_2)$ are true.
\subsubsection{When $\lambda <0$}
\begin{example}
Consider the nonlinear three point SBVP
\begin{eqnarray}
\label{Ex4.1.3}&& -y''(x)-\frac{\alpha}{x} y'(x)=1-2y^3,\\
\label{Ex4.2.3}&&~~y'(0)=0,~~~~~~ y(1)=\frac{1}{3} y\left(\frac{1}{4}\right).
\end{eqnarray}
\end{example}
Here $f(x,y)=1-2y^3,~~\delta=\frac{1}{3},~\eta =\frac{1}{4}$ and $\alpha$ satisfies $(H'_0)$. In this problem we choose lower and upper solutions as $v_0= -1$ and $u_0=1$, where $v_0\leq u_0$, i.e., it is well ordered case. The nonlinear term satisfies all assumptions for Theorem \ref{Theorem3.3}, and Lipschitz constant is $M_3=6$. Now we can find out a range for $\lambda<0$ such that the conditions $M_3+\lambda \leq 0$ and $(H'_0)$ are true.
\begin{example}
Consider the nonlinear three point SBVP
\begin{eqnarray}
\label{Ex4.1.4}&& -y''(x)-\frac{3}{x} y'(x)=1-7y^2,\\
\label{Ex4.2.4}&&~~y'(0)=0,~~~~~~ y(1)=2.2 y\left(\frac{1}{5}\right).
\end{eqnarray}
\end{example}
Here $f(x,y)=1-7y^3,~~\delta=2.2,~\eta =\frac{1}{5}$ and $\alpha$ satisfies $(H'_0)$. In this problem we choose lower and upper solutions as $v_0= 0$ and $u_0=\frac{7}{4}+\frac{5}{2}x^2$ i.e., it is well ordered case. The nonlinear term satisfies all assumptions for Theorem \ref{Theorem3.3}, and Lipschitz constant is $M_3=\frac{119}{2}$. Now we can find out a range for $\lambda<0$ such that the conditions $M_3+\lambda \leq 0$ and $(H'_0)$ are true.