Does a closed form solution to 1-D heat diffusion equation with Neumann and convective Boundary conditions exist? this problem has been eluding me and I've begun to wonder if a solution exists, and hoping I've simply overlooked something.
Using the normalized diffusion equation I have 
$$
{\partial{\rm T}\left(x,t\right) \over \partial t}
=
{\partial^{2}{\rm T}\left(x,t\right) \over \partial x^{2}}
$$
subjected to the boundary conditions:
$$
\partial_{x^*} T(0,t^*) = \dfrac{qL}{k(T(0,0)-T_\inf)}
$$
$$
\partial_{x^*} T(1,t^*) = -\dfrac{hL}{k}T^*(1,t^*)
$$
$$
u(x^*,0)=T_0
$$
where q,L,k,h are parameters to the system. Also, for all my intentions, $T_0$ is constant
Is this a well posed question? Over-constrained? Are there any known (semi-)analytical solutions?
I've attempted solutions using Laplace transform, which is non-invertible; Separation of variables, which leaves a residual constant (after solving for the eigenvalues) without an equation to solve for it; and integral approximation techniques, for which I have developed a quadratic polynomial solution. I'd be more than happy to share matlab code.
~~* edit ~~*
I'll try to explain my process a bit further, and exempand upon the main problem I'm having with the separation of variables approach-- satisfying the boundary conditions. If $T^*$ is broken into $X(x^*)\phi(t^*)$, then $\frac{X''}{X} = \frac{\dot{\phi}}{\phi} = -\lambda^2$ which gives $\phi = c_1*e^{(-\lambda^2*t)}$ and $X(x^*) = D_1sin(\lambda x^*) + D_2cos(\lambda x^*)$. From here, I can solve the boundary conditions to give:
$$\lambda = \frac{q L}{k D_1(T(0,0)-T_{inf})}$$
and $$D_1 = D_2Q \quad$$
$$ Q = \frac{\lambda sin(\lambda) - \frac{hL}{k}cos(\lambda)}{\lambda cos(\lambda) + \frac{hL}{k}sin(\lambda)} $$
So, $\lambda$ is transcendental and coupled to $D_2$. $D_2$ is (I think) indeterminable. I could just set $D_2=1$, but I don't think that's kosher. If this is the case, then the system is determinable, and $c_1$ can be found by orthogonality/Fourier series. If $D_2$ must be kept non-singular, then I'm at a loss.
Thanks so much!
-DH
 A: Okay, not giving up hope due to encouragement from @uvs I decided to snoop around the idea more. It turns out that in order to solve two inhomogeneous Neumann BCs, one must use Green's functions. Using this library from University of Nebraska, Lincoln, I was able to work out the solution for this post. I've had to change the problem to be
$$\dfrac{1}{\alpha}{\partial{\rm G}\left(x,t\right) \over \partial t}
=
{\partial^{2}{\rm G}\left(x,t\right) \over \partial x^{2}}$$
$$\partial_{x^*} G(0,t) = 0$$
$$\partial_{x^*} G(1,t) = 0$$
with the initial condition
$$T(x,0)=T_0(x)$$
and the forcing function
$$F(x') =p\space  \delta(x')\space H(t')  $$
where p is the power inserted upon the tissue ( which is equal to $q\space A/k$), $\delta$ is the dirac delta function, and H is the heaviside distribution function, and A is the x-sectional area of the material. This corresponds to problem type X23 in the GF library, and is solvable by the convolution
$$ T(x,t) = \int_0^t{\int_0^L{G(x,t|x',t')F(x',t') dx'}dt'} $$
which provides me an answer too ugly to be represented here, but I've solved this in matlab and arrived at some good results: here are some examples of the medium with 
${
h = 900;
k = 2.5;
rho = 998;
L = 0.0024;
Cp = 860;
alpha = k/(rho*Cp);
Tb = 36;
Tinf = 20;
T0 = Tb;
tf = 10;
D = 0.01;
Area = .25*pi*(D*D);
q = 5;
q = q/Area;
Bi = h*L/k;
absorptionCoeff = 0.5$


A: $\newcommand{\+}{^{\dagger}}%
 \newcommand{\angles}[1]{\left\langle #1 \right\rangle}%
 \newcommand{\braces}[1]{\left\lbrace #1 \right\rbrace}%
 \newcommand{\bracks}[1]{\left\lbrack #1 \right\rbrack}%
 \newcommand{\ceil}[1]{\,\left\lceil #1 \right\rceil\,}%
 \newcommand{\dd}{{\rm d}}%
 \newcommand{\down}{\downarrow}%
 \newcommand{\ds}[1]{\displaystyle{#1}}%
 \newcommand{\equalby}[1]{{#1 \atop {= \atop \vphantom{\huge A}}}}%
 \newcommand{\expo}[1]{\,{\rm e}^{#1}\,}%
 \newcommand{\fermi}{\,{\rm f}}%
 \newcommand{\floor}[1]{\,\left\lfloor #1 \right\rfloor\,}%
 \newcommand{\half}{{1 \over 2}}%
 \newcommand{\ic}{{\rm i}}%
 \newcommand{\iff}{\Longleftrightarrow}
 \newcommand{\imp}{\Longrightarrow}%
 \newcommand{\isdiv}{\,\left.\right\vert\,}%
 \newcommand{\ket}[1]{\left\vert #1\right\rangle}%
 \newcommand{\ol}[1]{\overline{#1}}%
 \newcommand{\pars}[1]{\left( #1 \right)}%
 \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
 \newcommand{\pp}{{\cal P}}%
 \newcommand{\root}[2][]{\,\sqrt[#1]{\,#2\,}\,}%
 \newcommand{\sech}{\,{\rm sech}}%
 \newcommand{\sgn}{\,{\rm sgn}}%
 \newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}}
 \newcommand{\ul}[1]{\underline{#1}}%
 \newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert}$
Use the Laplace transform $\ds{\tilde{\rm T}\pars{x,s} = \int_{0}^{\infty}{\rm T}\pars{x,t}\expo{-st}\,\dd t}$:

$$
\pars{\partiald[2]{}{x} - s}\tilde{\rm T}\pars{x,s} = - T_{0}
$$

$$
\tilde{\rm T}\pars{x,s} =
{T_{0} \over s} + A\sinh\pars{\root{s}\bracks{x - \half}}
+ B\cosh\pars{\root{s}\bracks{x - \half}} 
$$
Determine $A$ and $B$ with the Laplace transforms of $\partial_{x}{\rm T}\pars{x,t}$ conditions at $x = 0$ and $x = 1$.

With $\tilde{\rm T}\pars{x,s}$ you can determine ${\rm T}\pars{x,t}$. 

