First off, I'd note that your initial condition doesn't satisfy the boundary conditions, so you might want to instead use $u_0(x) = e^{-x^2} - e^{-L^2}$.
A great sanity check for problems like yours is the conservation property -- the total mass of $u$ should stay the same.
$$\begin{align}
\frac{d}{dt}\int_{-L}^Lu\, dx & = \int_{-L}^L\frac{\partial u}{\partial t}dx \\
& = \int_{-L}^L\frac{\partial}{\partial x}\left(vu + u^2\frac{\partial u}{\partial x}\right)dx \\
& = \left(vu + u^2\frac{\partial u}{\partial x}\right)\Big|_{x=-L}^{x=L} \\
& = 0,
\end{align}$$
because you've assumed that $u = 0$ at both endpoints.
Here I've written $v = \text{sign}(x)$ for the advection field, but this relation would hold true regardless of what $v$ was as long as you had the same nonlinear diffusion coefficient.
From the plot you've shown, it looks as if the numerical solutions you're obtaining are monotonically decreasing, which would violate the conservation property.
That suggests that there's an error in your numerical implementation somewhere.
When I run into problems like these, I usually try and come up with a simpler system and see if I can solve that first.
For example, what happens if you take out the advection term?
The PDE
$$\partial_tu = \partial_x(u^2\partial_xu)$$
is challenging enough by itself -- it's a free boundary problem.
Similarly, what happens if you take out the diffusion term and then smooth over the advection field?
Can you get a good approximation to the solutions of
$$\partial_tu = \partial_x(\tanh(x/\epsilon)u)$$
for different values of $\epsilon$?
Start with $\epsilon = L / 2$ and then see how things go as you decrease it to be equal to the mesh spacing $\delta x$.
You might even be able to write down an analytical solution using the method of characteristics.
Both of these simplified problems have conservation principles and other intrinsic mathematical properties that you can use as sanity checks.
