This kind of problem can be managed using a Fourier series for $x \in (0,1)$ and writing the general solution
$$
u(x,t)=\sum_{n=-\infty}^{+\infty}c_n(t)e^{i2\pi nx}
$$
and $u(x,0)=\sum_{n=-\infty}^{+\infty}c_n(0)e^{i2\pi nx}$. So, consider the equation
$$
\partial_tu(x,t)=D(x)\partial_{xx}u(x,t)
$$
and insert the given solution into the equation. You will get a set of ODEs for the coefficients $c_n(t)$ as
$$
\dot c_n(t)=\sum_{m}c_m(t)(-4\pi^2m^2)D_{nm}
$$
where I have put
$$
D_{nm} = \int_0^1dxD(x)e^{i2\pi(m-n)x}.
$$
The problem is then reduced to a proof of existence and uniqueness of a set of ODEs making it somewhat more manageable and reducing to conditions for the matrix of the diffusion coefficient.
Proof of uniqueness: Let us consider a generic matrix $D_{mn}$, we write $c_n(t)=c_n(0)e^{\lambda_nt}$ and substitute into the equations for the coefficients $c_n$. We get the eigenvalue problem
$$
(D-\lambda I)\bar c=0.
$$
being $D$ the matrix formed by the elements $-4\pi m^2D_{mn}$ and $\bar c$ the vector formed by the $c_n$ coefficients. Then we can state the following:
Theorem: If the given eigenvalue problem has a solution with all distinct eigenvalues, this is unique given the initial condition $u(x,0)=f(x)$ that fixes the coefficients $c_n(0)$ unequivocally.
This is true, provided the function $f(x)$ admits a Fourier series expansion $f(x)=\sum_nc_n(0)e^{i2\pi nx}$.
This theorem does not give an existence proof. This can be obtained by perturbation theory. So, let us turn back to the set of ODEs we started from. We can rewrite them as
$$
\dot c_n(t)=c_n(t)\epsilon_n+\sum_{m\ne n}c_m(t)(-4\pi^2m^2)D_{nm}
$$
where I have set $\epsilon_n=-4\pi^2n^2D_{nn}$. We can rewrite this as
$$
\dot b_n(t)=\sum_{m\ne n}b_m(t)e^{(\epsilon_m-\epsilon_n)t}(-4\pi^2m^2)D_{nm}.
$$
Now, one has $u(x,0)=f(x)=\sum_nc_n(0)e^{i2\pi nx}$ and using $c_n(0)=b_n(0)$ one can iterate the above set of equations to get
$$
b_n(t)=c_n(0)+\sum_{m\ne n}(-4\pi^2m^2)D_{nm}c_m(0)\int_0^te^{(\epsilon_m-\epsilon_n)t'}dt'
$$
$$
+\sum_{m\ne n}(-4\pi^2m^2)D_{nm}\sum_{k\ne m}(-4\pi^2k^2)D_{mk}c_k(0)\int_0^tdt'e^{(\epsilon_m-\epsilon_n)t'}\int_0^{t'}dt''e^{(\epsilon_k-\epsilon_m)t''}+\ldots
$$
It is worth nothing the following: 1) If this series converges then the solution exists and is unique in agreement with the theorem provided above. 2) For $t\rightarrow\infty$ this series contains polynomial terms in $t$, starting from second order, that can harm it. The point 2) shows that secular terms can be present and need to be resummed: they just give higher order corrections to the eigenvalues giving a series for them that must converge to the eigenvalue $\lambda_n$ introduced in the uniqueness theorem.
When $D$ is just a constant, we recover from all this the ordinary solution of this parabolic equation as it should and all converges nicely. This means that, for any $D(x)$ a check for convergence must be performed using the above series and also that degeneracies are not present in the eigenvalue problem.