Diffusion with spatially variable coefficient and source I've carefully studied many resources (such as 1, 2, 3, 4 and etc.) where they discuss different forms of diffusion equations. However, my specific equation was not studied in there. I'm hoping to get help on solving a diffusion equation with spatially variable coefficient and source, as follows (1D): 
$$\frac{\partial c}{\partial t}=\frac{\partial}{\partial x}J=\frac{\partial}{\partial x}(D\frac{\partial c}{\partial x} +E)$$
$$\sigma(x,0)=\sigma_0 \\
J(-L/2,t)=J_{-} \\
J(+L/2,t)=J_{+} $$
where $-L/2\leq x \leq L/2$ and, $D(x)$ and $E(x)$ are only spatially variable.
Separation of variables and Laplace's transforms are both welcome. Approximate form of solution with series analysis are highly appreciated.

More details as requested:
The original equation to solve is as follows ($-L/2\leq x \leq L/2$):
$$\frac{\partial{\sigma}}{\partial t}=\frac{\partial}{\partial x}J=B\frac{\partial}{\partial x}(-\frac{D}{kT}eZ\frac{\partial{V}}{\partial x}+\frac{D}{kT}\frac{Q}{T}\frac{\partial{T}}{\partial x}-\frac{D}{kT}\Omega\frac{\partial{\sigma}}{\partial x})$$
$$\sigma(x,0)=\sigma_0 \\
J(-L/2,t)=J_{-} \\
J(+L/2,t)=J_{+} $$
where 


*

*Diffusivity: $D(x)=D_0e^{-\frac{E_a}{kT}}$ ☹️

*Voltage: $V(x)=-\rho j x$ 

*Temperature: $T(x)=T_0+T_m \big( 1-\cosh(\frac{L}{2\Gamma})\cosh(\frac{x}{\Gamma}) \big) +T_n \big( \sinh(\frac{L}{2\Gamma})\sinh(\frac{x}{\Gamma}) \big)$ 

*Mean hydro-static stress: $\sigma(x,t)$: 


and constants are


*

*$B=3e10[Pa]$: Bulk modulus

*$k=1.38e-23[J/K]$: Boltzmann constant

*$e=1.6e-19[C]$: Electron charge

*$Z=1$: Charge number

*$Q=1.6e-19[J]$: Specific heat transfer

*$\Omega=1.66e-29[m^3]$: Atomic volume

*$D_0=7.56e-5[m^2/s]$ : Diffusion prefactor

*$E_a=1.6e-19[J]$ : Activation energy

*$\rho=1.67e-8[ohm.m]$ : Resistivity

*$T_0=\frac{T_-+T_+}{2}$ : Approx external temperature 

*$T_n=\frac{T_--T_+}{2}$ : Ends Temperature diff 

*$T_m=\frac{j^2\rho \Gamma^2}{k_{th}}$ : Max Temperature rise

*$\Gamma =6e-6[m] $ : Heating characteristic length

*$k_{th}=400[W/m.K]$ : Thermal conductance


But really, you only need to know the following inputs (roughly center ranges 
 added):


*

*$L\sim 100[um]$ : Wire length

*$j\sim 5e10[A/m^2]$ : Current density

*$T_-=T(x=-L/2)\sim373[K]$ : Fixed temperature at $-L/2$ end 

*$T_+=T(x=+L/2)\sim373[K]$ : ixed temperature at $+L/2$ end

*$\sigma(x,0) = \sigma_0$


The original equation can be seen as:
$$\sigma_t=(C_1\frac{D}{T}+C_2\frac{D}{T^2}T_x+C_3\frac{D}{T}\sigma_x)_x$$
By further encapsulation, it can be seen as follows (similar to what mentioned in the question):
$$\sigma_t=(F\sigma_x)_x+E_x$$
 A: I would suggest the introduction of suitable functions $u$ and $v$ such that $c = u + v$ and $u$ satisfies the non-homogeneous boundary conditions
\begin{align*}
x = -L/2: & \quad u_x = (J_o - E_o)/D_o, \\
x =  L/2: & \quad u_x = (J_L - E_L)/D_L,
\end{align*} where the subscripts for $D$, $J$ and $E$ indicate evaluation at the corresponding endpoint. The easiest function $u$ is $u(x) = \frac{1}{2} A ( x - B)^2$, where the constants $A$ and $B$ may be determined from the two conditions above.
Now, notice that the problem for $v$ is homogeneous in the spatial directions and given by
$$ v_t = ( D v_x)_x + Q(x), \quad -L/2 < x < L/2, $$
with BCs $v_x(-L/2,t) = v_x(L/2,t) = 0$ and IC $v(x,0) = c(x,0) - u(x)$. The source term is $Q(x) = (D u_x)_x + E_x$, which is also known. 
Thanks to superposition, your problem has now homogeneous boundary conditions but contains a nonzero, in general, source term. You may find now nontrivial solutions for the problem with $Q = 0$ (eigenfunctions), and expand the solution to the problem in terms of these functions (basis). 
Can you take it from here?

Follow up:
Let's consider for a moment the problem $w_t = (D w_x)_x$ with homogenous boundary conditions $w_x(-L/2,t) = w_x(L/2,t) = 0$. I don't care too much about the initial conditions for now. Since $D = D(x)$, the problem is separable and admits a combination of solutions of the form $w = X(x)T(t)$, where $X$ and $T$ are nonzero. Introduction of the ansatz yields
$$ X T' = T \, (D X')' \implies \frac{T'}{T} = \frac{1}{X} (D X')' = \lambda, $$ where the primes denote differentiation with respect the corresponding independent variable and $\lambda$ is a real (negative, positive or zero) constant. Now we need some information about $D$, since the shape of the final solution will strongly depend on the definition of $D$. Is it positive everywhere in the domain? Is it linear/quadratic/cubic...? 
Note that the problem in $x$ is 
$$ \frac{\mathrm{d}}{\mathrm{d}x} \left( D \frac{\mathrm{d}X}{\mathrm{d
 }x} \right) - \lambda X = 0, \quad -L/2 < x < L/2, \quad X'(-L/2) = X'(L/2) = 0, $$ where $\lambda$ is a priori unknown (it has to be chosen appropriately as an eigenvalue). This problem may not have an easy to obtain closed-form solution. For $D = 1$ (or any other positive constant) the solution take form of cosines and sines but for slightly more complicated choices of $D$ you arrive at Bessel, Airy, Legendre functions, etc. 
If you give me a hint on what kind of $D$s you're trying to solve the problem for, I can elaborate on finding a basis for your eigenfunctions on your domain and expand $v$ in terms of this basis (the $w$ functions).

Hope this helps.
