This can be solved numerically using implicit Euler method. The partial differential equation:
$$
\frac{\partial u}{\partial t}=D(x)\frac{\partial^2 u}{\partial x^2}+\frac{d D(x)}{d x}\frac{\partial u}{\partial x}
$$
Subjected to the following initial and boundary conditions:
\begin{align*}
u(x,0)&=\sin(\pi x)\\
u(0,t)&=0\\
u(1,t)&=0
\end{align*}
With finite difference, we can write:
$$
\frac{u_{i,j+1}-u_{i,j}}{\Delta t}=D(x_i)\big[\frac{u_{i+1,j+1}-2u_{i,j+1}+u_{i-1,j+1}}{\Delta x^2}\big]+\frac{d D(x)}{d x}\Big|_{x=x_i}\cdot\big[\frac{u_{i+1,j+1}-u_{i,j+1}}{\Delta x}\big]
$$
Note that $i$ refers to the $i$th node of the discretized space, and $j$ corresponds to the $j$th node of the discretized time. The equation can be rearranged, such that:
$$
u_{i,j+1}-D(x_i)\big[\frac{u_{i+1,j+1}-2u_{i,j+1}+u_{i-1,j+1}}{\Delta x^2}\big]\Delta t-\frac{d D(x)}{d x}\Big|_{x=x_i}\cdot\big[\frac{u_{i+1,j+1}-u_{i,j+1}}{\Delta x}\big]\Delta t=u_{i,j}
$$
Now let us write the terms in the following form:
\begin{align*}
A_{i,0}&=-D(x_i)\frac{\Delta t}{\Delta x^2}\\
A_{i,1}&=1+2D(x_i)\frac{\Delta t}{\Delta x^2}+\frac{d D(x)}{d x}\Big|_{x=x_i}\frac{\Delta t}{\Delta x}\\
A_{i,2}&=-D(x_i)\frac{\Delta t}{\Delta x^2}-\frac{d D(x)}{d x}\Big|_{x=x_i}\frac{\Delta t}{\Delta x}
\end{align*}
Therefore, at each time step we solve this equation:
$$
\begin{bmatrix}
A_{1,1} & A_{1,2} & 0 & 0&\cdots & 0 & 0 & 0\\
A_{2,0} & A_{2,1} & A_{2,2} & 0&\cdots & 0 & 0 & 0\\
0 & A_{3,0} & A_{3,1} & A_{3,2}&\cdots & 0 & 0 & 0\\
\vdots & \vdots & \vdots & \ddots & \ddots & \ddots & \vdots& \vdots \\
0 & 0 & 0 & 0&\cdots& A_{N-3,0} & A_{N-3,1} & A_{N-3,2} \\
0 & 0 & 0 & 0&\cdots & 0 & A_{N-2,0} & A_{N-2,1}
\end{bmatrix}\begin{bmatrix}
u_{1,j+1}\\u_{2,j+1}\\u_{3,j+1}\\ \vdots \\u_{N-3,j+1}\\u_{N-2,j+1}
\end{bmatrix}=\begin{bmatrix}
u_{1,j}\\u_{2,j}\\u_{3,j}\\ \vdots \\u_{N-3,j}\\u_{N-2,j}
\end{bmatrix}
$$
With $u_0$ and $u_{N-1}$ always equal to zero. Assume $D(x)$ is a sigmoidal function:
$$
D(x)=\frac{1}{1+e^{-x}}
$$
$$
\frac{dD(x)}{dx}=\frac{e^{-x}}{(1+e^{-x})^2}
$$
We then implement this using Python 3.5. The program is as follows:
import numpy as np
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings("ignore")
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.cm as cm
def f1(x):
return 1/(1+np.exp(-x))
def f2(x):
return np.exp(-x)/(1+np.exp(-x))**2
dx=0.01
x=np.linspace(0.01,0.99,99)
x2=np.linspace(0,1,101)
dt=0.001
t_step=500
sum_mat=np.zeros((int(t_step/25+1),len(x)+2))
u=np.zeros(len(x))
for i in range(len(x)):
u[i]=np.sin(np.pi*x[i])
t=[]
for j in range(t_step):
if j==0:
t.append(j*dt)
c_count=0
for i in range(len(x)):
sum_mat[c_count][i+1]=u[i]
c_count=c_count+1
else:
Am=np.zeros((len(x),len(x)))
for i in range(len(x)):
if i==0:
Am[i][i]=1+f1(x[i])*2*dt/dx**2+f2(x[i])*dt/dx
Am[i][i+1]=-f1(x[i])*dt/dx**2-f2(x[i])*dt/dx
elif i==len(x)-1:
Am[i][i]=1+f1(x[i])*2*dt/dx**2+f2(x[i])*dt/dx
Am[i][i-1]=-f1(x[i])*dt/dx**2
else:
Am[i][i]=1+f1(x[i])*2*dt/dx**2+f2(x[i])*dt/dx
Am[i][i-1]=-f1(x[i])*dt/dx**2
Am[i][i+1]=-f1(x[i])*dt/dx**2-f2(x[i])*dt/dx
u_next=np.linalg.solve(Am,u)
u=u_next
if (j+1)%25==0:
t.append(j*dt)
for i in range(len(x)):
sum_mat[c_count][i+1]=u_next[i]
u[i]=u_next[i]
c_count=c_count+1
X,T=np.meshgrid(x2,t)
fig = plt.figure()
plt.rc('font', **{'family': 'serif', 'serif': ['Computer Modern']})
plt.rc('text', usetex=True)
plt.gcf().set_size_inches(4,3,forward=True)
plt.subplots_adjust(left=0.18,bottom=0.15)
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X,T,sum_mat,cmap=cm.jet,lw=0.1,rstride=1, cstride=1)
plt.xlabel("$x$")
plt.ylabel("$t$")
ax.set_zlabel('$u$')
fig.savefig('ans.png',transparent=True,dpi=300)
Finally, this link is the plot of the numerical solution.
