Since $c$ doesn't play any role, denote $u_\infty(x)=e^{-V(x)}$, also $D=(0,1)\times(0,\infty)$ is the domain in which the solution is sought. Put $p(x,t)=u(x,t)/u_\infty(x,t)=e^{V(x)}u(x,t)$. The problem for the new function $p$ becomes
$$
\begin{cases}
\partial_t p =e^{V(x)} \partial_x \left(u_\infty \partial_x p\right)= \partial_x^2p+V'(x)\partial_xp&\text{ in } D \\
\partial_x p= 0 &\text{ for } t>0, x \in \{ 0,1 \} \\
p|_{t=0} =e^{V(x)} u_0(x) &\text{ in } [0,1]
\end{cases}
$$
The equation is uniformly parabolic in $D$ with smooth enough coefficients and the corresponding results apply.
Nothing is said about the initial function $u_0$ so I will assume for now it belongs to $C^{2+\alpha}([0,1])$, where $\alpha\in(0,1)$. If additionally to require the compatibility conditions (between the initial and boundary data) to be fulfilled, in this case
$$
(e^{V(x)}u_0(x))'|_{x=0}= (e^{V(x)}u_0(x))'|_{x=1}=0,
$$
then there exists a unique solution $p$ which belongs to the anisotropic Holder space $C^{2+\alpha,1+\alpha/2}(\bar D)$, see, for example O. A. Ladyzhenskaya, V.A. Solonnikov, N.N. Uraltseva, Linear and quasi-linear equations of parabolic type, ch. 4.
Returning to the initial problem one gets
- $u\in C^{2,1}(\bar D)$.
Also I think that the solution is smooth wrt $t$ in $\bar D_\varepsilon=D\cap\{t>\varepsilon\}$, $\varepsilon>0$, namely,
$\partial_t^mu\in C(\bar D_\varepsilon)$ for any $m\in \mathbb N$, though I haven't a reference.
If to reduce the regularity of $u_0$, say $u_0\in L_1(0,1)$, then the assertions in $\bar D_\varepsilon$ still hold: $u\in C^{2,1}(\bar D_\varepsilon)$, $\partial_t^mu\in C(\bar D_\varepsilon)$.
- If $u_0\in C([0,1])$ then $u$ is continuous up to the moment $t=0$: $u\in C(\bar D)$.
- If $u_0$ is non-negative then $u$ is non-negative.
- $\lim_{t\to+\infty}u(x,t)=Ke^{-V(x)}$, where
$$
K=\left(\int_0^1u_\infty\,dx\right)^{-1} \int_0^1u_0\,dx.
$$
Rationale for pp. 5,6 is as following. The boundary conditions of the second kind for $p$ mean physically that the endpoints are isolated and heat doesn't flow through them. In the absence of heat sources the initial distribution $u_0$ is flattened out to the mean "temperature" $K$, i.e. $\lim_{t\to+\infty}p(x,t)=K$ and $\lim_{t\to+\infty}u(x,t)=Ke^{-V(x)}$. Note that it is a stationary solution of the original equation.
To find $K$ we'll integrate the equation and use the boundary conditions:
$$
\int_0^t\int_0^1\partial_tu(x,\tau)\,dxd\tau=
\int_0^1u(x,t)\,dx-\int_0^1u(x,0)\,dx=
$$
$$
\int_0^t\int_0^1
\partial_x \left(u_\infty \partial_x \left( \frac{u}{u_\infty} \right)\right)\,dxd\tau=
\int_0^tu_\infty \partial_x \left( \frac{u}{u_\infty} \right)\Bigg|_{x=0}^{x=1}\,d\tau=0,
$$
so
$$
\int_0^1u(x,t)\,dx=\int_0^1u(x,0)\,dx.
$$
Taking the limit $t\to+\infty$ gives
$$
\int_0^1K e^{-V(x)}\,dx=\int_0^1u_0(x)\,dx.
$$
Update
Rewriting the equation in the form
$$
\partial_tp=u_\infty^{-1}\partial_x(u_\infty\partial_xp),
$$
the problem for $p$ can be solved by the Fourier method (separation of variables) using the Sturm–Liouville theory. Namely, denote
$$
(f,g)=\int_0^1u_\infty fg\,dx
$$
a scalar product with weight $u_\infty$ and $\|f\|=(f,f)^{1/2}$.
Let $\{\psi_n\}$, s.t. $\|\psi_n\|=1$, and $0=\lambda_1<\lambda_2<\ldots$ be an orthonormal system of eigenfunctions and eigenvalues of the corresponding Sturm–Liouville problem.
The first eigenfunction (for $\lambda_1=0$) is constant: $\psi_1=\|u_\infty\|^{-1/2}$.
Then for $u_0\in L_2(0,1)$
$$
p(x,t)=\sum_{n=1}^\infty a_n\psi_n(x)e^{-\lambda_n t}=K+\sum_{n=2}^\infty a_n\psi_n(x)e^{-\lambda_n t},
$$
where $a_n=\frac{(\psi_n,u_0)}{\|\psi_n\|}$.
From there we have p.6. and also p.2 because the series for $\partial_t^mp$ converge for $t>0$. It follow from this representation that the convergence rate of $u$ to $Ku_\infty$ is exponential (at least as $e^{-\lambda_2t}$).
