I have this non-linear partial differential equation. $$ \frac{\partial C}{\partial t}=\left(\frac{\partial C}{\partial x}\right)^2+C\frac{\partial^2 C}{\partial x^2} $$
I want to use the finite difference method to solve it either with the implicit method or the Crank-Nicolson method, witch I have done with linear PDE's, but how is this done with non-linear PDE's
First off, the PDE can be rewritten instead as
$$\frac{\partial C}{\partial t} = \frac{\partial}{\partial x}C\frac{\partial C}{\partial x}$$
or, by applying the product rule in reverse again, as
$$\frac{\partial C}{\partial t} = \frac{1}{2}\frac{\partial^2}{\partial x^2}C^2.$$
This equation is often referred to as the porous medium equation or the slow diffusion equation and it's got some fascinating properties. The regular old linear heat equation
$$\frac{\partial C}{\partial t} = \frac{\partial}{\partial x}D\frac{\partial C}{\partial x}$$
with conductivity $D$ has the property that disturbances propagate with infinite speed. The slow diffusion equation, on the other hand, has a well-defined front that propagates with finite speed, which is very unusual. This one of the simplest examples of a free boundary problem.
There are many ways to discretize nonlinear partial differential equations. I'll describe a very blunt and simple approach that should hopefully get you started; it's based on using the first form of the slow diffusion equation that I wrote down. The idea is to pretend like it's the linear diffusion equation, but to then take the diffusion coefficient to be the value of $C$ from the previous timestep. Let $$C_k^n = C(k\cdot\delta x, n\cdot\delta t)$$ be the values of the solution at time step $n$ and and spatial point $k$ and let $D_{k+\frac{1}{2}}$ be the value of the diffusion coefficient in the interval $[k\cdot\delta x, (k + 1)\delta x]$. An implicit Euler discretization of the linear diffusion equation would be
$$\frac{C^{n + 1}_k - C^n_k}{\delta t} = \frac{D_{k + \frac{1}{2}}\frac{C^{n + 1}_{k + 1} - C^{n + 1}_k}{\delta x} - D_{k - \frac{1}{2}}\frac{C^{n + 1}_k - C^{n + 1}_{k - 1}}{\delta x}}{\delta x}.$$
You can then take the diffusion coefficient in each interval as
$$D_{k + \frac{1}{2}} = \frac{C_{k + 1}^n + C_k^n}{2}$$
using the concentration from the previous timestep to approximate the nonlinearity.
If you want a more accurate numerical solver, you might want to look into implementing Newton's method. That will be more involved though and it's better to try simple approaches first.
An important thing to account for is that parabolic PDE have a maximum principle -- a solution that is initially positive will stay positive. The implicit Euler discretization has its own discrete maximum principle, whereas the Crank-Nicholson method does not. Despite the fact that the Crank-Nicholson method is formally of higher order in $\delta t$, it could actually be much worse in practice for this reason.
