If I understand the question I think you should in fact be able to use an implicit method such as the Crank-Nicolson. Whenever such a method is used, we do not know the value at the next time-step, but we set up a system of equations which we can solve to find the value. The non-linearity you have does not prevent you from using an implicit method but it might make it slow to do so, as you will possibly have to invert an $N$ by $N$ matrix at each timestep (Where $N$ is the number of lattice points in your spatial discretisation).

---

In your FP example you have some operator which acts on $u$ to give you the time derivative. Let's give that operator a symbol to make things easier to understand:

$$ F(u) equiv frac{partial^2}{partial x^2}left(D(b(u(x,t)))u(x,t)right) $$

The Crank-Nicolson method is then:

$$frac{u^{t+1} - u^t}{dt} = frac{1}{2}left(F^{t+1}(u) + F^{t}right(u))$$

Where the superscripts denote which the timestep (and I haven't written in anything to do with the spatial position but it will be involved in the finite difference scheme for the spatial derivatives). I take it your confusion arises from the term $F^{t+1}(u)$? You do not know it at this timestep so how can you use it? The answer is that you don't straightforwardly use it, it is part of what you are calculating. By setting up the finite difference scheme you will end up with a matrix with loads of unknown elements in it, but this is fine. You take all of those elements over to the left hand side and then solve for everything all at once using some linear algebra software. This is what makes the method implicit, the fact that you have to invert this matrix.

---

Edit:

In order to solve for the time evolution of $D$ it depends on how exactly $D$ behaves.

Option 1

Your original post made it sound as if it is simply a function of $u$, (via some variable $b(u)$), perhaps something like $D = au(t)^2$ or something. If this is the case then simply solve for $u(t)$ and you can calculate $D$.

Option 2

If however $D$ has its own distinct time dynamics then they will need to be specified. You didn't mention that you had an expression for $dD/dt$ so I assumed you were in the option 1 case. If you have an expression for the time evolution then you can then set up a CN scheme for $D$ too - combining this with the one you have for $u$ will give you a $2N$ by $2N$ matrix to invert. You will end up with something vaguely like:

$$begin{pmatrix} P & Q R & S end{pmatrix} begin{pmatrix} u^{t+1} D^{t+1} end{pmatrix} = begin{pmatrix} G(u^{t}) H(D^{t}) end{pmatrix}$$

Although the exact details of the matrices involved will depend on the form of the time evolution.

Option 3

If you do not know exactly $D$ varies in time and it can't be calculated from other variables then you cannot know how the system will evolve.

You can use an implicit scheme if you iterate each time step.

• Start by estimating $u^{t+1}$ somehow from the previously computed data - at worst, you could take $u^{t+1} = u^t$, or (probably better, after you have done the first time step) extrapolate from $u^t$ and $u^{t-1}$.

• Use the estimated value of $u^{t+1}$ to calculate an estimate for $D^{t+1}$.

• Do a Crank-Nicholson step using the estimate for to $D^{t+1}$ to get an updated estimate for $u^{t+1}$.

• Repeat until $u^{t+1}$ converges.

I don't have any experience with the Fokker-Planck equation, but for a mildly nonlinear equation and/or small enough time steps the iteration will often converge in 2 or 3 passes. Since the time steps for an implicit method are usually much bigger than those for an explicit method, to get the same accuracy, doing a few iterations per step is still a win in terms of computation time.

The next level of sophistication is to take variable-sized time steps, depending on how the solution is behaving. If the general form of the solution is some sort of "transient" behaviour depending on the boundary conditions, superimposed on a "steady state" solution, you can use what you know about the "transient" behaviour to help to choose a good way to vary the time steps. For example in the classic application of Crank-Nicholson to a convection-diffusion equation, the "transient" part of the solution is typically a decaying exponential function of time, and successive time steps can increase in a geometric progression without losing accuracy.

For severely nonlinear problems, it may be worth using this iterative approach with a backward-difference scheme in time, which will improve the numerical stability of the solution at the expense of reducing the accuracy. This can be an improvement over a forward-difference (explicit) solution in situations where the forward-difference time steps are limited by stability, not by accuracy.