\documentclass[12p]{article}
\usepackage{graphicx}
\usepackage{amssymb,amsmath}
\usepackage[pdftex,a4paper,bookmarks=true,hyperindex=true,
colorlinks=true]{hyperref}
\usepackage[final]{pdfpages}
\usepackage{amssymb,dsfont}
\newtheorem{Theorem}{Theorem}
\newtheorem{Remark}{Remark}
\newtheorem{Lemma}{Lemma}
\newtheorem{Corollary}{Corollary}
\newtheorem{Definition}{Definition}
\numberwithin{Definition}{section} \numberwithin{Theorem}{section}
\numberwithin{Lemma}{section}
\newtheorem{Example}{Example}
%\numberwithin{equation}{section}
\numberwithin{Corollary}{section} \numberwithin{Example}{section}
\numberwithin{Remark}{section}
\usepackage{epstopdf}
\begin{document} 
\title{ Implicit finite difference approximation for time fractional heat conduction under boundary condition of second kind } 
\author{ \small T.N. Mishra*, \small K.N. Rai
\\ \small \textit{DST-CIMS, Faculty of Science, BHU, Varanasi}\\\textit{*Corresponding author E-mail: {t.mishra01@gmail.com}} }
\date{}
\maketitle 
\begin{abstract}
The time fractional heat conduction in an infinite plate of finite thickness, when both faces are subjected to boundary conditions of second kind, has been studied. The time fractional heat conduction equation is used, when attempting to describe transport process with long memory, where the rate of heat conduction is inconsistent with the classical Brownian motion. The stability and convergence of this numerical scheme has been discussed and observed that the solution is unconditionally stable. The whole analysis is presented in dimensionless form. A numerical example of particular interest has been studied and discussed in details.
\end{abstract}
\begin{small}
\textbf{ Keyword :} Caputo fractional derivative, Time fractional heat conduction, Boundary condition of second kind, Implicit finite difference scheme, Unconditional stability, Convergence.
\end{small}
\section{Introduction}
\par Fractional differential equations have attracted increasing attention because they have applications in various field of science and engineering. Many phenomena in engineering and applied sciences can be described successfully by developing model using fractional calculus, such as material and mechanics, signal processing,  anomalous diffusion, biological system, finance, hydrology and so on (cf. \cite{podlubny:book1999},\cite{meer:icassp2001},\cite{solomon:icassp1993},\cite{yuste:icassp2001},\cite{raberto:icassp2002},\cite{mainardi:icassp2000},\cite{ben:icassp2000}). The main physical purpose for investigating diffusion equations of fractional order is to describe phenomena of anomalous diffusion in transport processes through complex and/or disordered system including fractal media and fractional kinetic equations. Fractional kinetic equations have proved particularly useful in the context of anomalous slow diffusion (sub-diffusion).
\par Fractional sub-diffusion equation is a class of anomalous diffusive system, which is obtained by replacing the time derivative in ordinary diffusion by a fractional derivative of order $ \alpha $ with $ 0<\alpha<1 $ (in Riemann-Liouville or Caputo sense). It is a known fact that the anomalous diffusion is characterized by a diffusion constant and the mean square displacement of diffusing species in the form
\begin{equation*}
\langle x^{2}(t)\rangle\sim t^{\alpha}, \hspace{1.00 cm}t\rightarrow \infty\hspace{8.0 cm}
\end{equation*}
where, $\alpha$ $(0<\alpha<1)$ is the anomalous diffusion exponent. 
\par Numerical approaches to solve different types of fractional diffusion models have been increasingly appearing in literature. The stochastic solution of space-time fractional diffusion equation was proposed by Meerchaert \cite{meer:icassp2002}. Huang and Liu \cite{huang:icassp2005}, has obtained the solution of time fractional diffusion equation and advection-dispersion equation in $n$-dimensional whole space and in half space respectively, in terms of Green function by Schneider and Wyss. Yuste and Acedo \cite{yuste:icassp2005} and Yuste \cite{yuste:icassp2006}, presented an explicit scheme and weighted average finite difference scheme, respectively for the time fractional diffusion equation and analyzed the stability of these two schemes by von Neumann method. Zhuang and Liu \cite{zhuang:icassp2006}, approximated the time fractional diffusion equation by implicit difference scheme. Based on Grunwald-Letnikov formula, Chen and Liu \cite{chen:icassp2007}, constructed the difference scheme  for the fractional sub-diffusion equation, and showed the stability and convergency of the difference scheme using the Fourier method.
Zhuang et al. \cite{zhuang:icassp2008}, obtained an implicit numerical method for solving sub-diffusion equation by integrating the original equation on the both sides. Stability and convergency of the scheme were proved by the energy method. Murio \cite{murio:icassp2008}, presented implicit finite difference approximation for time fractional diffusion equations but in showing stability of the scheme by Fourier method he made some flaw. Ding and Zhang \cite{ding:icassp2011}, showed the stability of implicit finite difference approximation for time fractional diffusion equations.
Chen et al. \cite{chen:icassp2008}, constructed an implicit and an explicit finite difference methods for the fractional reaction-sub-diffusion equation and analyzed the stability and convergency of schemes by Fourier method. Zhang \cite{zhang:icassp2009}, constructed the unconditionally stable finite difference method for solving fractional partial differential equation. Singh et al.\cite{singh:icassp2011}, solved the fractional bioheat equations governing the process of heat transfer in tissues during thermal therapy by finite difference method and homotopy perturbation method.

\par The work mentioned above are dealing with Dirichlet boundary conditions, where no boundary discretization errors are involved. However, for Neumann boundary condition (i.e boundary condition of second kind), the discretization of boundary conditions must dealt with carefully to match global accuracy. Langlands and Henry \cite{lang:icassp2005}, developed an implicit difference scheme with convergence order $o(\tau+h^{2})$, based on $L_{1}$ approximation, for Riemann-Liouville fractional derivative and numerically verified the unconditional stability of difference scheme but without global convergence analysis. Recently, Zhao and Sun \cite{zhao:icassp2011}, proposed a box-type scheme for solving a class of fractional sub-diffusion equation with homogeneous Neumann boundary conditions.
Since many problems in science and engineering involve Neumann boundary conditions \cite{lang:icassp2005}, \cite{povst:icassp2010},\cite{naber:icassp2004}, 
such as zero flow or specified flow flux condition. Thus, it is very desirable to construct high-order algorithms for efficient computations. This motivates us to consider the implicit finite difference approximation for spatial discretization.
\par The main purpose of this paper is to construct a mathematical model for time fractional heat conduction in an infinite plate of finite thickness and discuss the anomalous diffusion, sub-diffusion and effect of Fourier number $(F_{0})$, Predvoditelev number $(P_{d})$ on heat conduction. Further, the asymptotic behavior of time fractional constant $\alpha$ on heat transport has to be studied.
The rest of the paper is organized as follows. Section $1$ deals fundamental of fractional calculus, time fractional heat conduction equation and mathematical model of heat conduction. An implicit finite difference scheme is given in section $2$. Stability and convergence of the scheme is given in section $3$. Section $4$ contains  results and discussion. The conclusion and future research plan is given in section $5$.\\\\
\textbf{Nomenclature}
\begin{tabbing}
\hspace{7.5cm}\=\kill
$a_{q}$\>thermal diffusivity $(m^{2} s^{-1})$ \\
$b_{1}$\> time dependent coefficient ($s^{-1}K)$\\
$ c_{b}$ \>specific heat capacity (J/kg  $K$)\\
$ _{a}D^{\alpha}_{t}  $ \>Grunwald Letnikove fractional derivative \\ 
$ _{a}^{RL}D^{\alpha}_{t}$ \>Riemann-Liouville fractional derivative \\
$ _{a}^{C}D^{\alpha}_{t}$\>Caputo fractional derivative \\
g(t)\> a continuous differentiable function of $t$  \\
$h_{1}$\>scale factor \\ 
$K$\>thermal conductivity ($W/m$ $K$)\\
$k$\>temporal grid size\\
$q(t)$\> heat flux ($W/$ $m^{2}$)  \\
$Q$\>heat source $(J)$ \\ 
2R\>thickness of the infinite plate $(m)$ \\
$T$\>temperature $(K)$\\
$T_{0}$\>initial temperature $(K)$ \\
$T_{c}$\>ambient temperature $(K)$ \\
$\triangledown T$\> temperature gradient ($K$/m) \\
$\Delta x$\> spatial grid size ($=h$)\\
$\alpha$\> order of fractional derivative in (n-1, n] \\ 
$\rho$ \> density (kg $m^{-3})$   
\end{tabbing}
\textbf{Dimensional variable and similarity criteria} 
\begin{tabbing}
\hspace{7.5 cm}\=\kill
$f(x)$\>dimensionless initial temperature \\ 
$F_{0} $\>Fourier number, $\frac{a_{q}t}{R^{2}}$\\
$k_{i}(F_{0})$ \>Kirpichev number,  $\frac{q(\tau)R}{K (T_{c}-T_{0})}$\\
$k_{i_{1}}(F_{0})$ \>dimensionless heat flux at boundary  x=-1 \\
$k_{i_{2}}(F_{0})$\>dimensionless heat flux at boundary x= 1\\
$P_{0}(x,F_{0})$\>Pomerantsev number, $\frac{Q(x,\tau)R^{2}}{K (T_{c}-T_{0})}$ \\
$P_{d}$ \>Predvoditelev number, $\frac{b_{1} R^{2}}{a_{q}\Delta T}$ \\
$x$ \>dimensionless spatial coordinate, $\frac{r}{R}$\\
$\theta$\>dimensionless temperature, $\frac{T-T_{0}}{T_{c}-T_{0}}$
\end{tabbing}
\subsection{Fundamental of fractional calculus}
The Grunwald Letnikov fractional derivative $ _{a}D^{\alpha}_{t} g(t)$ \cite{podlubny:book1999}, of order $ \alpha > 0$, of a function $ g(t) $ is defined as
\begin{equation}\label{1.1}
 _{a}D^{\alpha}_{t} g(t)=\sum\limits_{k=0}^{m}\frac{g^{k}(a)(t-a)^{(-\alpha+k)} }{\Gamma (-\alpha +k+1)}+\frac{1}{\Gamma (-\alpha +k+1)}\int\limits_{a}^{t}(t-\tau)^{m-\alpha} g^{m+1}(\tau) d\tau
\end{equation}
This definition of fractional derivative is defined under the assumption that functions are $ m+1 $ - times continuously differentiable in the closed interval [ a, t ] where, $ m $ is an integer number satisfying the condition $ m > \alpha-1 $. But such type of functions are very narrow.
\par Let [a, t] be a finite interval on the real axis $\Re$. The Riemann-Liouville (R-L) fractional derivative \cite{podlubny:book1999} $^{RL}_{a}D^{\alpha}_{t} g(t)  $, of order $\alpha > 0$, of a function $ g(t)$ is defined as
\begin{equation}\label{1.2}
 _{a}^{RL}D^{\alpha}_{t}g(t)=\frac{d^{n}}{d t^{n}}[\frac{1}{\Gamma (n-\alpha)}]\int\limits_{a}^{t}(t-\tau)^{n-\alpha-1} g(\tau) d\tau, \hspace{0.5cm} (n-1<\alpha\leq n)
\end{equation}
The Riemann-Liouville approach have played an important role in the development of theory of fractional calculus due to their applications in pure mathematics. However, R-L approach lead to initial conditions containing the limit values at the lower terminal $ t= a$ as
\begin{equation*}
\lim_{t \rightarrow a}\hspace{0.1 cm}  _{a}^{RL}D^{\alpha-j}_{t} g(t)=\bar {c_{j}}, \hspace{0.5cm} j=1,2,...n.
\end{equation*}
where, $\bar {c_{j}}$ are given constants. The above initial value problems with such initial conditions can be successfully solved mathematically
but their solutions have no known physical interpretation due to which they are useless.
\par
The Caputo fractional derivative $ _{a}^{C}D^{\alpha}_{t} g(t)$ \cite{podlubny:book1999} , of order $\alpha$, of a function $g(t)$ is defined as
\begin{equation}\label{1.3}
_{a}^{C}D^{\alpha}_{t} g(t)=\frac{1}{\Gamma (n-\alpha)}\int\limits_{a}^{t}(t-\tau)^{n-\alpha-1}g^{n}(\tau) d\tau,  \hspace{0.5cm} n-1<\alpha \leq n
\end{equation}
Under natural conditions on the function $g(t)$, for $\alpha \rightarrow n $ the Caputo derivative becomes a conventional n-th derivative of the function $g(t)$. The main advantage of Caputo's approach is that the initial conditions for fractional differential equations with Caputo derivatives takes on the same form as for integer-order differential equations, i.e. contain the limit values of integer-order derivatives of unknown functions at the lower terminal $t=a$, due to this in present study we take Caputo fractional derivative.
\subsection{Time fractional heat conduction equation }
The classical theory of heat conduction is based on Fourier law, 
$ \vec q = - \triangledown T $, relating to the heat flux vector $(\vec q)$ to the temperature gradient $(\triangledown T)$, where, thermal conductivity K is a geometry-independent
coefficient and mainly depends on the composition and structure of the material and the temperature. In many one dimensional systems with total momentum conservation, the heat conduction equation does not obey the Fourier law and heat conductivity depends on the system size \cite{khanna:book2004}. Such size dependent thermal conductivity has been observed in theoretical models, such as the harmonic chains \cite{rieder:icassp1967}, the FPU-β model \cite{lepri:icassp1998} and the hard point gas model \cite{grass:icassp2002}. \par The heat conduction equation is a combination of conservation equation of energy and conduction law. The energy conservation is expressed locally by
\begin{equation}\label{1.4}
-\triangledown\vec  q(t) = \rho c_{b}\frac{\partial T }{\partial t }
\end{equation}
where, $c_{b}$, is the specific heat capacity at constant pressure and  $\rho, $ is the density. The theory of heat conduction proposed by Norwood \cite{norw:icassp1972}, the corresponding generalization of Fourier law,
\begin{equation}\label{1.5}
\vec q(t)= -K \int\limits_{0}^{t}a (t-\tau)\triangledown T(\tau)d\tau
\end{equation}
where, $a (t-\tau)$ is the kernel.
The time non local dependence between the heat flux vector and the temperature gradient with long tale power kernel \cite{povst:icassp2008},
\begin{equation}\label{1.6}
\vec q(t)= \frac{-K}{\Gamma(\alpha)}\int\limits_{0}^{t} (t-\tau)^{\alpha-1}\frac{\partial \triangledown T (\tau)}{\partial t } d \tau ,\hspace{1cm} 0<\alpha \leqslant 1
\end{equation}
Considering the definition of Caputo fractional derivative for $0<\alpha\leq 1$
\begin{equation}\label{1.7}
_{0}^{c}D_{t}^{\alpha}g(t)=\frac{1}{\Gamma(1-\alpha)}\int\limits_{0}^{t}\frac { g^{1}(\tau)}{(t-\tau)^{\alpha}}d\tau,\hspace{1cm} 0<\alpha\leq 1
\end{equation}
Now Eq. (\ref{1.6}) becomes,
\begin{equation}\label{1.8}
\vec q(t)=-K(_{0}^{c}D_{t}^{\alpha-1}\triangledown T)
\end{equation}
Let the functional relation between orthogonal curvilinear coordinates $(u_{1},u_{2},u_{3})$ and rectangular coordinates $(x,y,z)$ be given as
\begin{equation*}
x=X(u_{1},u_{2},u_{3}),\hspace{0.7 cm}y=Y(u_{1},u_{2},u_{3}),\hspace{0.7 cm}z=Z(u_{1},u_{2},u_{3})
\end{equation*}
Then, the gradient of temperature in orthogonal coordinate system $(u_{1},u_{2},u_{3})$ \cite{jiang:icassp2010}, is given by 
\begin{equation*}
\bigtriangledown T=\bar{u_{1}}\frac{1}{h_{1}}\frac{\partial T}{\partial u_{1}}+\bar{u_{2}}\frac{1}{h_{2}}\frac{\partial T}{\partial u_{2}}+\bar{u_{3}}\frac{1}{h_{3}}\frac{\partial T}{\partial u_{3}}
\end{equation*}
Here, the coefficients $h_{1}$, $h_{2}$ and $h_{3}$ are called scale factor and can be evaluated by
\begin{equation}\label{1.9}
h^{2}_{i}=\left(\frac{\partial x}{\partial u_{i}}\right)^{2}+\left(\frac{\partial y}{\partial u_{i}}\right)^{2}+\left(\frac{\partial z}{\partial u_{i}}\right)^{2}, \hspace{0.3 cm} i=1,2,3
\end{equation} 
and $\bar{u_{1}}$, $\bar{u_{2}}$ and $\bar{u_{3}}$ be the unit direction vectors in the direction of $u_{1}$, $u_{2}$ and $u_{3}$ in the general orthogonal curvilinear coordinate system (for more about of orthogonal curvilinear coordinate see Appendix).
\par In present study, we take the one-dimensional rectangular coordinate system and by utilizing Eq. (\ref{1.9}), the scale factor $h_{1}$ is determined as
$h^{2}_{1}=1$ and hence the gradient of temperature in orthogonal coordinate system is given by
\begin{equation}\label{1.10}
\triangledown T =\bar{x}\frac{1}{h_{1}}\frac{\partial T}{\partial x} 
\end{equation}
Then, the heat flux vector $\vec{q}$ becomes,
\begin{equation}\label{1.11}
\vec{q}=-K(^{c}_{0}D^{1-\alpha}_{t}\bar{x}\frac{1}{h_{1}}\frac{\partial T}{\partial x})
\end{equation}
Therefore divergence of heat flux vector $\vec{q}$ in orthogonal curvilinear coordinate system is given by 
\begin{equation}\label{1.12}
\begin{split}
\triangledown \vec{q} = -\frac{1}{h_{1}}[\frac{\partial}{\partial x}K(_{0}^{c}D^{1-\alpha}_{t}\frac{1}{h_{1}}\frac{\partial T}{\partial x})] 
\end{split}
\end{equation} 
Then, the fractional heat conduction equation in orthogonal curvilinear coordinate system is given by substituting the result of Eq. (\ref{1.12}) into Eq. (\ref{1.4}) i.e.
\begin{equation*}
\frac{1}{h_{1}}[\frac{\partial}{\partial x}K(_{0}^{c}D^{1-\alpha}_{t}\frac{1}{h_{1}}\frac{\partial T}{\partial x})]=\rho c_{b}\frac{\partial T}{\partial t}
\end{equation*}
\subsection{Mathematical model of heat conduction in an infinite plate of finite thickness}
The transport process in an infinite plate of finite thickness $2R$, with the long memory where, the rate of heat conduction in the body is inconsistent with the classical Brownian motion model, is considered. At initial instant, the temperature of the body is a function of space coordinate and the faces  are subjected to boundary condition of second kind. The mathematical model describing this anomalous heat conduction in presence of internal heat source, in dimensionless form can be written as:
\begin{equation}\label{1.13}
\frac{\partial^{\alpha}\theta}{\partial F_{0}^{\alpha}} = \frac{\partial^{2}\theta}{\partial x^{2}}+P_{0}(x,F_{0})
\end{equation}
\begin{equation}\label{1.14}
\hspace{1.5 cm}(-1)^{j}\frac{\partial \theta (x, F_{0})}{\partial x}=k_{i_{j}}(F_{0}),\hspace{0.5 cm}    x=(-1)^{j} ,\hspace{0.2 cm} j=1,2
\end{equation}
\begin{equation}\label{1.15}
\theta(x,0)= f(x),\hspace{1cm} F_{0}>0
\end{equation}
\section{Implicit finite difference scheme}
To establish the numerical approximation scheme, let $\Delta x = h =\frac{2}{p} > 0 $ and $k= \frac{F_{0}}{M}$  be the grid size in space and time  direction respectively. The grid points in the space interval $[-1, 1]$ are the numbers $x_{i}=-1+ih ,i=0,1,2,...p $ and grid points in the time interval $ [0,  F_{0}]$ are labelled $ F_{0_{n}} = n k, n=0,1,2,..,M $. The values of the functions $\theta $ and $ f $ at grid points are denoted by $ \theta_{i}^{n}=\theta(x_{i},F_{0_{n}})$ and $ f_{i}=f(x_{i})$, respectively.
In the differential Eq. (\ref{1.13}), we have adopted a symmetric second difference quotient in space at level $F_{0}=F_{0_{n+1}}$, for approximating the second order space derivative. The time fractional derivative term can be approximated by the following scheme 
\begin{equation}
\frac{\partial^{\alpha}\theta(x_{i},F_{0_{n+1}})}{\partial F_{0}^{\alpha}}\approxeq\frac{1}{\Gamma(1-\alpha)}\int\limits_{0}^{F_{0_{n+1}}}\frac{\partial\theta(x_{i},s)}{\partial F_{0}}(F_{0_{n+1}}-s)^{-\alpha}ds \hspace{6.9 cm}
\nonumber
\end{equation}
\begin{equation}
~~~~~~~~~\approxeq\frac{1}{\Gamma(1-\alpha)}\sum\limits_{j=0}^{n}\int\limits_{jk}^{(j+1)k}[\frac{\theta^{j+1}_{i}-\theta^{j}_{i}}{k}]\left\lbrace (n+1)k-s\right\rbrace ^{-\alpha}ds
\nonumber
\end{equation}
\begin{equation}
\hspace{2.5cm}\approxeq\frac{1}{\Gamma(1-\alpha)}\frac{1}{(1-\alpha)}\sum\limits_{j=0}^{n}([\frac{\theta^{j+1}_{i}-\theta^{j}_{i}}{k}][(n-j+1)^{1-\alpha}-(n-j)^{1-\alpha}])k^{1-\alpha}\hspace{4.5cm}
\nonumber
\end{equation}
\begin{equation}
\hspace{2.5cm}\approxeq\frac{1}{\Gamma(1-\alpha)}\frac{1}{(1-\alpha)}\sum\limits_{j=0}^{n}([\frac{\theta^{n-j+1}_{i}-\theta^{n-j}_{i}}{k}][(j+1)^{1-\alpha}-(j)^{1-\alpha}])k^{1-\alpha}\hspace{4.5cm}
\nonumber
\end{equation}
\begin{equation}
\begin{split}
\approxeq\frac{1}{\Gamma(1-\alpha)}\frac{1}{(1-\alpha)}\frac{1}{k^{\alpha}}(\theta^{n+1}_{i}-\theta^{n}_{i})\hspace{9.8cm}\\\hspace{2.5cm}+\frac{1}{\Gamma(1-\alpha)}\frac{1}{(1-\alpha)}\frac{1}{k^{\alpha}}\sum\limits_{j=1}^{n}(\theta^{n-j+1}_{i}-\theta^{n-j}_{i})[(j+1)^{1-\alpha}-(j)^{1-\alpha}]\hspace{5.4 cm}
\end{split}
\nonumber
\end{equation}
By taking $\sigma_{\alpha,k}=\frac{1}{\Gamma(1-\alpha)}\frac{1}{(1-\alpha)}\frac{1}{k^{\alpha}}$ and shifting indices $w^{\alpha}_{j}=(j+1)^{1-\alpha}-(j)^{1-\alpha}$, and define
\begin{equation}
L^{\alpha}_{h,k}\theta(x_{i,F_{0_{n+1}}})=\sigma_{\alpha,k}\sum\limits_{j=0}^{n}w^{\alpha}_{j}(\theta^{n+1-j}_{i}-\theta^{n-j}_{i})
\nonumber
\end{equation}
We have,
\begin{equation}
|\frac{\partial^{\alpha}\theta(x_{i},F_{0_{n+1}})}{\partial F_{0}^{\alpha}}-L^{\alpha}_{h,k}\theta(x_{i,F_{0_{n+1}}})|\hspace{8.5 cm}
\nonumber
\end{equation}
\begin{equation}
\hspace{1.0 cm}\leqslant\frac{1}{(1-\alpha)}k\sum\limits_{j=0}^{n}\int\limits_{F_{0_{j}}}^{F_{0_{j+1}}}|\frac{\partial\theta(x_{i},\xi)}{\partial \xi}-\frac{\theta(x_{i},F_{0_{j+1}})-\theta(x_{i},F_{0_{j}})}{k}|\frac{d \xi}{(F_{0_{n+1}}-\xi)^{\alpha}}
\nonumber
\end{equation}
\begin{equation}
\leqslant\frac{C_{1}}{\Gamma(1-\alpha)}k \sum\limits_{j=0}^{n}\int\limits_{F_{0_{j}}}^{F_{0_{j+1}}}\frac{d \xi}{(F_{0_{n+1}}-\xi)^{\alpha}}
\nonumber
\end{equation}
\begin{equation}
\leqslant\frac{C_{1}}{\Gamma(1-\alpha)}k \int\limits_{0}^{F_{0_{j+1}}}\frac{d \xi}{(F_{0_{n+1}}-\xi)^{\alpha}}~~
\nonumber
\end{equation}
\begin{equation}\label{2.1}
|\frac{\partial^{\alpha}\theta(x_{i},F_{0_{n+1}})}{\partial F_{0}^{\alpha}}-L^{\alpha}_{h,k}\theta(x_{i,F_{0_{n+1}}})|\leqslant \bar C_{1}k\hspace{4.0cm}
\end{equation}
where, $C_{1} $ and $ \bar C_{1}$ are constants.
\par We now approximating derivatives at $x_{0}$ and $x_{p}$, by Newtons five point formula, Eq. (\ref{1.14}) becomes
\begin{equation}\label{2.3}
\theta_{0}=\frac{20hk_{i_{1}}+13\theta_{1}+17\theta_{2}-9\theta_{3}}{21},\hspace{0.2 cm}j=1 \hspace{1.2 cm}
\end{equation}
\begin{equation}\label{2.4}
\theta_{p}=\frac{20hk_{i_{2}}+13\theta_{p-1}+17\theta_{p-2}-9\theta_{p-3}}{21},\hspace{0.2 cm}j=2
\end{equation}
On using finite central difference formula and from Eqs. (\ref{2.3})-(\ref{2.4}), Eqs. (\ref{1.13})-(\ref{1.15}), reduces in the vector-matrix form as follows:
\begin{equation}\label{2.5}
A\Theta^{1}=\bar{\Theta}^{0}+ P^{1}\hspace{7.80cm}
\end{equation}
\begin{equation}\label{2.6}
A\Theta^{n+1}=c_{1}\bar{\Theta}^{n}+c_{2}\bar{\Theta}^{n-1}+............+c_{n}\bar{\Theta}^{1}+d_{n}\bar{\Theta}^{0}+ P^{n}, n\geqslant 1\hspace{1.8cm}
\end{equation}
\begin{equation}\label{2.7}
\Theta^{0}=f \hspace{8.70cm}
\end{equation}
where,\\
\begin{equation*}
\bar{\Theta}^{n}=\Theta^{n}\left(\sigma_{\alpha,k}\right)_{1\times 1}\hspace{10.5 cm}
\end{equation*}
\begin{equation*}
\Theta^{n}=\begin{pmatrix}
\theta^{n}_{1}&\theta^{n}_{2}&\theta^{n}_{3}&.&.&.&.&\theta^{n}_{p-2}&\theta^{n}_{p-1} 
\end{pmatrix}^{t}\hspace{6.3cm}
\end{equation*}
\begin{equation*}
 A=\begin{pmatrix}
 a_{1}&-38 r_{1}&9r_{1}&0&.&.&.&0 \\ -r&b&-r&0&.&.&.&0 \\ 0&-r&b&-r&.&.&.&0 \\ \\0&.&.&.&-r&b&-r&0\\0 &.&.&.&0&-r&b&-r\\0&.&.&.&0&9r_{1}&-38 r_{1}&a_{1}
\end{pmatrix}_{(p-1)\times(p-1)}\hspace{1.0cm}
\end{equation*}
\begin{equation*}
P_{0}^{n}=\begin{pmatrix}
P_{0_{1}}^{n}+20 h k_{i,1}r_{1}&P_{0_{2}}^{n}&P_{0_{3}}^{n}&.&.&.&.&P_{0_{p-3}}^{n}&P_{0_{p-2}}^{n} &P_{0_{p-1}}^{n}+20 h k_{i,2}r_{1}
\end{pmatrix}^{t}\hspace{4.8cm}
\end{equation*}
\begin{equation*}
f=\begin{pmatrix}
f(x_{1})&f(x_{2})&f(x_{3})&.&.&.&.&f(x_{p-3})&f(x_{p-2}) &f(x_{p-1})
\end{pmatrix}^{t}\hspace{3.0cm}
\end{equation*}
$c_{j}=(w^{\alpha}_{j-1}-w^{\alpha}_{j})$,\hspace{1.4 cm}  $d_{n}=w^{\alpha}_{n}$,\hspace{1.2 cm}
$a_{1}=(\sigma_{\alpha,k}+29 r_{1})$\\~~$b=(\sigma_{\alpha,k}+2 r)$,\hspace{2.0 cm}$r_{1}=\frac{1}{21 h^{2}}$,\hspace{1.3 cm} $r=\frac{1}{ h^{2}}$
\\here, $X^{t}$ means the transpose of matrix $ X $.\\
From Eqs. (\ref{2.5})-(\ref{2.7}) we observe that at each time step, we get a triangular system of linear equations which utilize all the history of the computed solution up to that time. In present analysis, Mathematical software MATLAB 7.0 has been used to obtain dimensionless temperature $\Theta^{n}$. 
\section{Stability and convergence}
In this section we will analyze the stability and convergence of the implicit finite difference scheme for time fractional heat conduction equation. From Eqs. (\ref{2.5})-(\ref{2.7}), and by using the property of the function $g(x)=x^{1-\alpha},(1\leqslant x)$ and
\begin{equation*}
1=w^{\alpha}_{0}>w^{\alpha}_{1}>w^{\alpha}_{2}>w^{\alpha}_{3}>w^{\alpha}_{4}>..... \longrightarrow 0
\end{equation*}
we have the following,
\begin{equation*}
\sum\limits_{j=1}^{n}c_{j}=1- w^{\alpha}_{n}\hspace{5.0 cm}
\end{equation*}
\begin{equation*}
\sum\limits_{j=1}^{\infty}c_{j}=1, 2-2^{1-\alpha}>2^{2-\alpha}-3^{1-\alpha}-1=c_{1}>c_{2}>c_{3} ....\longrightarrow 0
\end{equation*}
Hence, we have
\begin{Remark}
\begin{itemize}
\item 
In Eqs. (\ref{2.5})-(\ref{2.7}), Matrix A is strictly diagonally dominant with positive diagonal term and non positive off-diagonal terms. Hence the Eqs. (\ref{2.5})-(\ref{2.7}) can be solved.
\item The solution $\Theta^{n}$ (n=1,2,3....M) of Eqs. (\ref{2.5})-(\ref{2.7}), preserve non negativity, if $\Theta^{0}$ is non negative.
\end{itemize}
\end{Remark}
Now we suppose $\tilde{\theta^{j}_{i}}$ is the approximate solution of Eqs. (\ref{2.5})-(\ref{2.7}), then the error $\epsilon_{i}^{j}=\tilde{\theta^{j}_{i}}-\theta^{j}_{i}$ satisfies,
\begin{equation*}
\begin{split}
AE^{1}= \bar{E^{0}}\hspace{9.2cm}\\
AE^{n+1}=c_{1} \bar{E}^{n}+c_{2}\bar{E}^{n-1}+c_{3}\bar{ E}^{n-2}+............+c_{n} \bar{E}^{1}+d_{n}\bar{ E}^{0}, n\geqslant 1\hspace{0.7 cm}
\end{split}
\end{equation*}
where, $\bar{E}^{n}=E^{n}\left(\sigma_{\alpha,k}\right)_{1\times 1}$\\and 
~~~~$E^{n}=\begin{pmatrix}
\epsilon^{n}_{1}&\epsilon^{n}_{2}&\epsilon^{n}_{3}&.&.&.&.&\epsilon^{n}_{p-2}&\epsilon^{n}_{p-1}
\end{pmatrix}^{t}$.\\ Hence, by the mathematical induction following result can be proved.
\begin{Theorem}
 $||E^{n}||_{\infty}\leqslant ||\bar{E}^{0}||_{\infty}$, n=1,2,3....
\end{Theorem} 
\textbf{Proof :} For n=1,
\begin{equation*}
-r\epsilon_{i-1}^{1}+(\sigma_{\alpha,k}+2 r)\epsilon_{i}^{1}-r \epsilon_{i+1}^{1}=\sigma_{\alpha,k}\epsilon_{i}^{0}
\end{equation*}
Let $|\epsilon_{l}^{1}|=\underbrace{max}_{1\leqslant i \leqslant p-1}|\epsilon_{i}^{1}|$ then we have
\begin{equation*}
|\epsilon_{l}^{1}|\leqslant -r|\epsilon_{l-1}^{1}|+(\sigma_{\alpha,k}+2 r)|\epsilon_{l}^{1}|-r|\epsilon_{l+1}^{1}|
\end{equation*}
\begin{equation*}
~~~\leqslant|-r\epsilon_{l-1}^{1}+(\sigma_{\alpha,k}+2 r)\epsilon_{l}^{1}-r \epsilon_{l+1}^{1}|
\end{equation*}
\begin{equation*}
=|\sigma_{\alpha,k}\epsilon_{l}^{0}|\hspace{3.5 cm}
\end{equation*}
\begin{equation*}
~~||E^{1}||_{\infty}\leqslant||\bar{E}^{0}||_{\infty}\hspace{5.0 cm}
\end{equation*}
Therefore, ~~~~~~$||E^{1}||_{\infty}\leqslant||\bar{E}^{0}||_{\infty}$.\\ 
Suppose that~ $||E^{j}||_{\infty}\leqslant||\bar{E}^{0}||_{\infty}$,\hspace{0.2 cm} j= 1,2,3....n. Let $|\epsilon_{l}^{n+1}|=\underbrace{max}_{1\leqslant i \leqslant p-1}|\epsilon_{i}^{n+1}|$, 
then we also have,
\begin{equation*}
|\epsilon_{l}^{n+1}|\leqslant -r|\epsilon_{l-1}^{n+1}|+(\sigma_{\alpha,k}+2 r)|\epsilon_{l}^{n+1}|-r|\epsilon_{l+1}^{n+1}|
\end{equation*}
\begin{equation*}
~~~~~\leqslant|-r\epsilon_{l-1}^{n+1}+(\sigma_{\alpha,k}+2 r)\epsilon_{l}^{n+1}-r \epsilon_{l+1}^{n+1}|
\end{equation*}
\begin{equation*}
=|c_{1}\epsilon_{l}^{n}+\sum\limits_{j=1}^{n-1}c_{j+1}\epsilon_{l}^{n-j}+d_{n}\epsilon_{l}^{0}|~~
\end{equation*}
\begin{equation*}
~~\leqslant c_{1}|\epsilon_{l}^{n}|+\sum\limits_{j=1}^{n-1}c_{j+1}|\epsilon_{l}^{n-j}|+d_{n}|\epsilon_{l}^{0}|
\end{equation*}
\begin{equation*}
||E^{n+1}||_{\infty}\leqslant c_{1}||E^{n}||_{\infty}+\sum\limits_{j=1}^{n-1}c_{j+1}||E^{n-j} ||+w^{\alpha}_{n}||E^{0}||_{\infty}
\end{equation*}
\begin{equation*}
\leqslant\left\lbrace c_{1}+\sum\limits_{j=1}^{n-1}c_{j+1} +w^{\alpha}_{n}\right\rbrace ||\bar{E}^{0}||_{\infty}
\end{equation*}
\begin{equation*}
=||\bar{E}^{0}||_{\infty}\hspace{4.5 cm}
\end{equation*}
\begin{equation*}
||E^{n+1}||_{\infty}\leqslant||\bar{E}^{0}||_{\infty}\hspace{6.5 cm}
\end{equation*}
Hence, the following,
\begin{Remark}
The time fractional implicit difference approximation defined by (\ref{2.5}) and (\ref{2.7}) are unconditionally stable.
\end{Remark}
\par Now we will discuss the convergence of implicit finite difference approximation for time fractional heat conduction equation under non-homogeneous boundary conditions of second kind. Let $\theta(x_{i},F_{0_{n}})$ be the exact solution of the time fractional heat conduction Eqs. (\ref{1.13})-(\ref{1.15}) at mesh point $(x_{i},F_{0_{n}})$. Define $e_{i}^{n}=\theta(x_{i},F_{0_{n}})-\theta^{n}_{i}, i=1,2,3.....p-1; n=1,2,3......M $ and $e^{n}=(\sigma_{\alpha,k}e^{n}_{1},\sigma_{\alpha,k}e^{n}_{2},\sigma_{\alpha,k}e^{n}_{3}.....\sigma_{\alpha,k}e^{n}_{p-1})^{t}$. Using
$e^{0}=0$, and substitution into (\ref{2.5})-(\ref{2.7}) leads to,
\begin{equation*}
-re_{i-1}^{1}+(\sigma_{\alpha,k}+2 r)e_{i}^{1}-re_{i+1}^{1}=R_{i}^{1}
\end{equation*}
\begin{equation*}
-re_{i-1}^{n+1}+(\sigma_{\alpha,k}+2 r)e_{i}^{n+1}-re_{i+1}^{n+1}=c_{1} e_{i}^{n}+\sum\limits_{j=1}^{n-1}c_{j+1} e_{i}^{n-j}+R_{i}^{n+1}
\end{equation*}
where,
\begin{equation*}
\begin{split}
R_{i}^{n+1}=\sigma_{\alpha,k}\theta(x_{i},F_{0_{n+1}})-\sigma_{\alpha,k}\theta(x_{i},F_{0_{n}})+\sigma_{\alpha,k}\sum\limits_{j=1}^{n-1}w^{\alpha}_{j}\left[ \theta(x_{i},F_{0_{n+1-j}})\right.\\\left.-\theta(x_{i},F_{0_{n-j}})\right]-r\left[ \theta(x_{i+1},F_{0_{n+1}})-2\theta(x_{i},F_{0_{n+1}})+\theta(x_{i-1},F_{0_{n+1}})\right]-P_{0_{i}}^{n+1}
\end{split}
\end{equation*}
\begin{equation}\label{3.1}
\begin{split}
=\sigma_{\alpha,k}\sum\limits_{j=0}^{n}w^{\alpha}_{j}\left[ \theta(x_{i},F_{0_{n+1-j}})-\theta(x_{i},F_{0_{n-j}})\right]-r\left[ \theta(x_{i+1},F_{0_{n+1}})\right.\\\left.-2\theta(x_{i},F_{0_{n+1}})+\theta(x_{i-1},F_{0_{n+1}})\right]-P_{0_{i}}^{n+1}\hspace{3.5 cm}
\end{split}
\end{equation}
From equation (\ref{2.1}), we have
\begin{equation}\label{3.2}
\hspace{1.0 cm}\sigma_{\alpha,k}\sum\limits_{j=0}^{n}w^{\alpha}_{j}\left[ \theta(x_{i},F_{0_{n+1-j}})-\theta(x_{i},F_{0_{n-j}})\right]=\frac{\partial^{\alpha}\theta(x_{i},F_{0_{n+1}})}{\partial F^{\alpha}_{0_{n+1}} }+\bar C_{1}k
\end{equation}
and 
\begin{equation}\label{3.3}
\frac{\theta(x_{i+1},F_{0_{n+1}})-2\theta(x_{i},F_{0_{n+1}})+\theta(x_{i-1},F_{0_{n+1}})}{h^{2}}=\frac{\partial^{2}\theta(x_{i},F_{0_{n+1}})}{\partial x^{2} }+C_{2}h^{2}~~
\end{equation}
Therefore from Eqs. (\ref{3.1}), (\ref{3.2}) and (\ref{3.3}) we get,
\begin{equation*}
|R_{i}^{n+1}|\leqslant C(k^{1+\alpha}+k^{\alpha}h^{2}), i=1,2,3...p-1; n=0,1,2,3...M.
\end{equation*}
where, C is a constant.
\begin{Theorem}
$||e^{n}||_{\infty}\leqslant C (w^{\alpha}_{n-1})^{-1}(k^{1+\alpha}+k^{\alpha}h^{2}), n=1,2,....M$ where, $||e^{n}||_{\infty}= \underbrace{max}_{1\leqslant i \leqslant p-1} |e_{i}^{n}|$ and C is a constant.
\end{Theorem}
\textbf{Proof :} We will prove the above theorem by mathematical induction. For $n=1$, let $||e^{1}||_{\infty}=|\sigma_{\alpha,k}e_{l}^{1}|=|\sigma_{\alpha,k}||e_{l}^{1}|=\underbrace{max}_{1\leqslant i \leqslant p-1}|\sigma_{\alpha,k}e_{i}^{1}|=|\sigma_{\alpha,k}|\underbrace{max}_{1\leqslant i \leqslant p-1}|e_{i}^{1}| $ i.e $||e^{1}||_{\infty}=|e_{l}^{1}|=\underbrace{max}_{1\leqslant i \leqslant p-1}|e_{i}^{1}|$, hence we have
\begin{eqnarray*}
~~~~~|e_{l}^{1}|\leqslant -r|e_{l+1}^{1}|+(\sigma_{\alpha,k}+2 r)|e_{l}^{1}|-r|e_{l-1}^{1}|\\\leqslant|-r e_{l+1}^{1} +(\sigma_{\alpha,k}+2 r) e_{l}^{1} -r e_{l-1}^{1} |~~\\=|R^{1}_{i}|\hspace{4.7 cm}\\\leqslant C (w^{\alpha}_{0})^{-1}(k^{1+\alpha}+k^{\alpha}h^{2})\hspace{1.8 cm}
\end{eqnarray*}
Suppose that $||e^{j-1}||_{\infty}\leqslant C (w^{\alpha}_{j})^{-1}(k^{1+\alpha}+k^{\alpha}h^{2})$, j=1,2,..n and $|\sigma_{\alpha,k}e_{l}^{n+1}|=\underbrace{max}_{1\leqslant i \leqslant p-1}|\sigma_{\alpha,k}e_{i}^{1}| $ i.e $|e_{l}^{n+1}|=\underbrace{max}_{1\leqslant i \leqslant p-1}|e_{i}^{1}|$.
Here we note that $ (w^{\alpha}_{j})^{-1}\leqslant (w^{\alpha}_{n})^{-1}$, j=0,1,2,3.....n. Then we have
\begin{equation*}
|e_{l}^{n+1}|\leqslant -r|e_{l+1}^{n+1}|+(\sigma_{\alpha,k}+2 r)|e_{l}^{n+1}|-r|e_{l-1}^{n+1}|
\end{equation*}
\begin{equation*}
\hspace{0.8 cm}\leqslant|-r e_{l+1}^{n+1} +(\sigma_{\alpha,k}+2 r) e_{l}^{n+1} -r e_{l-1}^{n+1}|
\end{equation*}
\begin{equation*}
=|c_{1}e_{l}^{n}+ \sum\limits_{j=1}^{n-1}c_{j+1}e_{l}^{n-j}+R^{n+1}_{i}|\hspace{0.4 cm}
\end{equation*}
\begin{equation*}
\hspace{1.5 cm}\leqslant c_{1}|e_{l}^{n}|+ \sum\limits_{j=1}^{n-1}c_{j+1}|e_{l}^{n-j}|+C (k^{1+\alpha}+k^{\alpha}h^{2})
\end{equation*}
\begin{equation*}
\hspace{2.6 cm}\leqslant c_{1}||e^{n}||_{\infty}+ \sum\limits_{j=1}^{n-1}c_{j+1}||e^{n-j}||_{\infty}+C (k^{1+\alpha}+k^{\alpha}h^{2})
\end{equation*}
\begin{equation*}
\hspace{2.2 cm}\leqslant\left[c_{1}+\sum\limits_{j=1}^{n-1}c_{j+1}+w^{\alpha}_{n}\right](w^{\alpha}_{n})^{-1}C (k^{1+\alpha}+k^{\alpha}h^{2})
\end{equation*}
\begin{equation*}
=(w^{\alpha}_{n})^{-1}C (k^{1+\alpha}+k^{\alpha}h^{2})\hspace{1.2 cm}
\end{equation*}
Since,
\begin{equation*}
\hspace{0.3 cm}\lim_{n \rightarrow \infty}\frac{(w^{\alpha}_{n})^{-1}}{n^{\alpha}}= \lim_{n \rightarrow \infty}\frac{n^{-\alpha}}{(n+1)^{1-\alpha}-n^{1-\alpha}}
\end{equation*}
\begin{equation*}
\hspace{1.7 cm}= \lim_{n \rightarrow \infty}\frac{n^{-1}}{(1+\frac{1}{n})^{1-\alpha}-1}
\end{equation*}
\begin{equation*}
\hspace{1.0 cm}= \lim_{n \rightarrow \infty}\frac {n^{-1}} {(1-\alpha)n^{-1}}
\end{equation*}
\begin{equation*}
~~~=\frac{1}{1-\alpha}\hspace{1.2 cm}
\end{equation*}
Hence, there is a constant $\bar{C}$,such that\\
$||e^{n}||_{\infty}\leqslant \bar{C} n^{\alpha}(k^{1+\alpha}+k^{\alpha}h^{2})$.\\ Since, $nk\leqslant F_{0_{n}}$ is finite, we obtain the following result.
\begin{Theorem}
Let $\theta^{n}_{i}$ be the approximate value of $\theta(x_{i},F_{0_{n}})$ computed by use of the difference scheme (\ref{2.5}-(\ref{2.7}) then there is a positive constant C such that\\ $|\theta^{n}_{i}-\theta(x_{i},F_{0_{n}})|\leqslant C (k+h^{2}), i=1,2,3...p-1, n=1,2,3...M.$
\end{Theorem}
\section{Results and discussion}
\par This section presents complete analysis of anomalous diffusion, sub-diffusion and also studies the effect of Fourier number $(F_{0})$ and Predvoditelev number $(P_{d})$ on temperature and heat flux. Further, the effect of limiting case $\alpha \rightarrow 0$ on temperature and heat flux has studied. The figures presented in this study, only the parameters whose values different from the reference value are indicated. The selected reference values include $h =0.2$, $k = \frac{1}{64}$ and $k_{i_{1}}(F_{0})=k_{i_{2}}(F_{0})=\exp(-P_{d}F_{0})$, $P_{0}=0$.\\\\
\includegraphics[width=6.5 cm,height=6.5 cm]{F:/MATLAB1/Figure3/F1a}
\includegraphics[width=6.5 cm,height=6.5 cm]{F:/MATLAB1/Figure3/F1b}\\
\begin{small}
\textbf{Fig. 1.} Variation of dimensionless temperature ($\Theta$) with space ($x$)
\end{small}\\\par
Variation of dimensionless temperature $\Theta$ with space coordinate $x$ is given in Fig. 1. The temperature decreases as $x$ increases and attains a minimum value at centre i.e $x=0$ and then increases symmetrically, which shows the parabolic behavior of heat conduction.\\\\
\includegraphics[width=6.5 cm,height=6.5 cm]{F:/MATLAB1/Figure3/F2a}
\includegraphics[width=6.5 cm,height=6.5 cm]{F:/MATLAB1/Figure3/F2b}\\
\begin{small}
\textbf{Fig. 2.} Variation of dimensionless temperature ($\Theta$) with $\alpha$
\end{small}
\\\par From Fig. 2, it is evident that, dimensionless temperature increases and attains a maximum value with increase in time fractional order $\alpha$.\\\\
\includegraphics[width=6.5 cm,height=6.0 cm]{F:/MATLAB1/Figure3/F3a}
\includegraphics[width=6.5 cm,height=6.0 cm]{F:/MATLAB1/Figure3/F3b}\\
\begin{small}
\textbf{Fig. 3.} Variation of dimensionless temperature ($\Theta$) with $\alpha$
\end{small}\\\par
A close examination of Figs. 2 and 3 reveal that,
\begin{itemize}
\item For $F_{0}\leqslant 0.375 $ (say critical $F^{*}_{0}$) and $P_{d}=1.0$, temperature response exhibit fast diffusion for any small time fractional order $\alpha$ at $x=0$.
\item For $F_{0}>0.375$ and $P_{d}=1.0$, sub-diffusion phenomena becomes apparent, consequently dimensionless temperature increases and attains a maximum value for certain time fractional order (say critical $\alpha^{*}$) and decreases further with increase of $\alpha$ as shown in Fig. 3(a).
\end{itemize}
\begin{center}
\includegraphics[width=8.5 cm,height=6.0 cm]{F:/MATLAB1/Figure3/F4}
\end{center}
\begin{small}
\textbf{Fig. 4.} Variation of $\alpha^{*}$ with $F_{0}$
\end{small}\\ 
\begin{itemize}
\item From Fig. 3(b) it is observe that as $P_{d}$ decreases, $F^{*}_{0}$ increases because Predvoditelev number characterizes the temperature time change of the outer environment i.e. Predvoditelev number $P_{d}=\frac{b_{1}t}{F_{0}(T_{c}-T_{0}})$.
\end{itemize}
\par When $F_{0}$ increases from $F^{*}_{0}$ ($=0.375$), for fixed $P_{d}=1.0$, the value of $\alpha^{*}$ at which maximum temperature attains, increases, as shown in Fig. 4.\\\\
\includegraphics[width=6.5 cm,height=6.0 cm]{F:/MATLAB1/Figure3/F5a}
\includegraphics[width=6.5 cm,height=6.0 cm]{F:/MATLAB1/Figure3/F5b}\\
\begin{small}
\textbf{Fig. 5.} Variation of dimensionless temperature ($\Theta$) with $F_{0}$
\end{small}\\
\par Variation of dimensionless temperature with $F_{0}$ for different $\alpha$ is given in Figs. 5 and 6. For very small time fractional order, temperature profile exhibit fast diffusion at fixed $P_{d}=1.0$, as shown in Fig. 5(a). From Figs. 5(b) and 6(a) it is evident that sub-diffusion occurs consequently temperature increases and attains a maximum value for certain $F^{*}_{0}$ and then decreases with increase of $\alpha$, similar to the result of Fig. 2.
\\\\
\includegraphics[width=6.5 cm,height=6.0 cm]{F:/MATLAB1/Figure3/F6a}
\includegraphics[width=6.5 cm,height=6.0 cm]{F:/MATLAB1/Figure3/F6b}\\
\begin{small}
\textbf{Fig. 6.} Variation of dimensionless temperature ($\Theta$) with $F_{0}$
\end{small}\\
\par When $\alpha\longrightarrow 1$, sub-diffusion phenomena disappear and for $\alpha=1$, temperature profile exhibit normal diffusion, as shown in Fig. 6(b).\\\\
\includegraphics[width=6.5 cm,height=7.0 cm]{F:/MATLAB1/Figure3/F7a}
\includegraphics[width=6.5 cm,height=7.0 cm]{F:/MATLAB1/Figure3/F7b}\\
\begin{small}
\textbf{Fig. 7.} Variation of dimensionless temperature ($\Theta$) with $P_{d}$
\end{small}
\\\par
Figs. 7 and 8 present temperature history with $P_{d}$. From Fig. 7, it is evident that, $P_{d}$ increases results $F_{0}$ decreases, consequently fast diffusion occurs for very small $\alpha$ close to zero, similar to result of Fig. 5(a) and sub-diffusion occurs for $\alpha > \alpha^{*}$ as shown in Fig. 8.\\\\
\includegraphics[width=6.5 cm,height=7.0 cm]{F:/MATLAB1/Figure3/F8a}
\includegraphics[width=6.5 cm,height=7.0 cm]{F:/MATLAB1/Figure3/F8b}\\
\begin{small}
\textbf{Fig. 8.} Variation of dimensionless temperature ($\Theta$) with $P_{d}$
\end{small}
\\\par
For very small $\alpha$ in the range of ($10^{-1},10^{-2}$), results fast diffusion, as shown in Fig. 9(a) and for $\alpha$ in the range ($10^{-1},10^{-10}$) i.e. $\alpha \rightarrow 0$, fractional tempera
\includegraphics[width=6.5 cm,height=7.0 cm]{F:/MATLAB1/Figure3/F9a}
\includegraphics[width=6.5 cm,height=7.0 cm]{F:/MATLAB1/Figure3/F9b}\\
\begin{small}
\textbf{Fig. 9.} Variation of dimensionless temperature ($\Theta$) with $\alpha$
\end{small}
\\\\-ture profile become asymptotically stable, for any fixed $P_{d}$ and any fixed $F_{0}$, as shown in Fig. 9(b).
\par Fig. 10 presents the dimensionless heat flux history at boundaries i.e. at locations $x=-1$ and $x=1$. On both boundary, heat flux exponentially decreases with $F_{0}$ for fixed $P_{d}$ and is constant with $\alpha$, as shown in Fig. 10(b).\\\\
\includegraphics[width=6.5 cm,height=7.0 cm]{F:/MATLAB1/Figure3/F11a}
\includegraphics[width=6.5 cm,height=7.0 cm]{F:/MATLAB1/Figure3/F11b}
\begin{small}
\textbf{Fig. 10.} Variation of dimensionless heat flux with $F_{0}$
\end{small}\\\par
Fig. 11 shows the executing time of our method in 2.3 GH i3 processor with 2 GB ram and Windows 7 operating system. 
\begin{center}
\includegraphics[width=8.5 cm,height=6.5 cm]{F:/MATLAB1/Figure3/F10}
\end{center}
\begin{small}
\textbf{Fig. 11.} Variation of CPU time with the time grid size
\end{small}
\section{Conclusions}
\par A mathematical model describing time fractional heat conduction (in the sense of Caputo), in an infinite plate of finite thickness whose faces are subjected to non-homogeneous boundary condition of second kind has been analyzed using a fully implicit unconditionally stable and convergent finite difference scheme.
\par Our simulation shows that temperature response exhibits fast diffusion when $F_{0}\leqslant F^{*}_{0}$, for very small $\alpha$ and any fixed $P_{d}$. When $F_{0}>F^{*}_{0}$, two cases arises, one is fast diffusion gradually disappear as $\alpha \longrightarrow \alpha^{*} $, and another is sub-diffusion becomes apparent for $\alpha>\alpha^{*}$. Further, for any fixed $F_{0}$ and $P_{d}$, when $\alpha \longrightarrow 0$, temperature profile becomes asymptotically stable. On both boundary heat flux exponentially decreases with $F_{0}$ at fixed $P_{d}$ and is constant with $\alpha$, for very small $\alpha$.
\par This method can be readily applied to two and three dimensional time as well as space fractional heat conduction problem. \\\\
\begin{small}
\textbf{ Acknowledgement } The authors of this article express their sincere thanks to the respected reviewers and editors for their valuable suggestions
in the improvement of the article. The authors are also thankful to Prof. Umesh Singh, Coordinator, DST-CIMS, BHU for providing the necessary facilities
during the manuscript.
\end{small}
\\\\
\textbf{Appendix: The orthogonal curvilinear coordinate} Let the functional relation between orthogonal curvilinear coordinates $(u_{1},u_{2},u_{3})$ and rectangular coordinates $(x,y,z)$ be given as
\begin{equation}\label{0.1}
x=X(u_{1},u_{2},u_{3}),\hspace{0.7 cm}y=Y(u_{1},u_{2},u_{3}),\hspace{0.7 cm}z=Z(u_{1},u_{2},u_{3})
\end{equation}
Then, the differential lengths $dx$, $dy$ and $dz$ are obtained from Eq. (\ref{0.1}) by differentiation
\begin{equation}\label{0.2}
dx=\sum\limits_{i=1}^{3}\frac{\partial X}{\partial u_{i}}du_{i}
\end{equation} 
\begin{equation}\label{0.3}
dy=\sum\limits_{i=1}^{3}\frac{\partial Y}{\partial u_{i}}du_{i}
\end{equation} 
\begin{equation}\label{0.4}
dz=\sum\limits_{i=1}^{3}\frac{\partial Z}{\partial u_{i}}du_{i}
\end{equation} 
The differential length $(dS)$ in orthogonal curvilinear coordinate system $(u_{1},u_{2},u_{3})$ is given by
\begin{equation}\label{0.5}
(dS)^{2}=h^{2}_{1}(du_{1})^{2}+h^{2}_{2}(du_{2})^{2}+h^{2}_{3}(du_{3})^{2}
\end{equation} 
where,\\
\begin{equation}
h^{2}_{i}=\left(\frac{\partial x}{\partial u_{i}}\right)^{2}+\left(\frac{\partial y}{\partial u_{i}}\right)^{2}+\left(\frac{\partial z}{\partial u_{i}}\right)^{2}, \hspace{0.3 cm} i=1,2,3
\end{equation}
Here, the coefficients $h_{1}$, $h_{2}$ and $h_{3}$ are called scale factor. Now the gradient of temperature in orthogonal coordinate system $(u_{1},u_{2},u_{3})$ is given by
\begin{equation}
\bigtriangledown T=\bar{u_{1}}\frac{1}{h_{1}}\frac{\partial T}{\partial u_{1}}+\bar{u_{2}}\frac{1}{h_{2}}\frac{\partial T}{\partial u_{2}}+\bar{u_{3}}\frac{1}{h_{3}}\frac{\partial T}{\partial u_{3}}
\end{equation}
where, $\bar{u_{1}}$, $\bar{u_{2}}$ and $\bar{u_{3}}$ be the unit direction vectors in the direction of $u_{1}$, $u_{2}$ and $u_{3}$ in the general orthogonal curvilinear coordinate system. To know more about the orthogonal curvilinear coordinate system, see literature \cite{ozisik:book1993}.\\

\bibliographystyle{elsart-num-sort}
\bibliography{appn}
\end{document}