I am trying to discretize the equation $u_t = \Delta (u^2)$ without using finite elements. If I use a Crank-Nicolson scheme I end up: $$\frac{u^{n+1}_j - u^n_j}{\Delta t} = \frac{(u^2)^{n+1}_{j+1} - 2(u^2)^{n+1}_j+(u^2)^{n+1}_{j-1}}{2(\Delta x)^2} + \frac{(u^2)^{n}_{j+1} - 2(u^2)^{n}_j+(u^2)^{n}_{j-1}}{2(\Delta x)^2} $$
Which then leads to $$u^{n+1}_j - \lambda ((u^2)^{n+1}_{j+1} - 2(u^2)^{n+1}_j+(u^2)^{n+1}_{j-1}) = u^n_j + \lambda((u^2)^{n}_{j+1} - 2(u^2)^{n}_j+(u^2)^{n}_{j-1}) $$ where $\lambda = \frac{\Delta t}{2 (\Delta x)^2}$.
I have no idea how to deal with the $(u^2)^{n+1}$ terms in my MATLAB code. Would it be best to use finite elements or is there a way that I can use finite difference methods for this problem?