I am trying to solve a very simple case of the chemical Fokker-Planch equations. It is basically a diffusion equation of the probability density function of the state of a chemical reaction network. The process is very simple (a single molecule, X, is produced at constant rate λ and degraded linearly at rate β).

I am just trying to solve it once, and create a plot of the time development. The process is very simple, so performance is not an issue (as long as it doesn’t take 2 months to solve…). However, I do want a correct(ish) solution.

The diffusion is in 1d space, and the equation looks like this:

∂P/∂t = - (∂/∂x)⋅[(λ-βx)⋅P] + (1/2)⋅(∂²/∂x²)⋅[(λ+βx)⋅P]

I tried to make a discrete version:

```
grid = 0:0.05:200.0; Δx = grid[2]
A(i,u,λ,β,grid) = (0 < i <= length(grid)) ? u[i]*(λ-β*grid[i]) : 0.0
B(i,u,λ,β,grid) = (0 < i <= length(grid)) ? u[i]*(λ+β*grid[i]) : 0.0
function cfpe(du,u,p,t)
λ,β = p
for i = 1:length(grid)
du[i] = u[i] - (0.5/Δx) * (A(i+1,u,λ,β,grid) - A(i-1,u,λ,β,grid)) + (0.5/Δx^2) * (B(i+1,u,λ,β,grid) -2*B(i,u,λ,β,grid) + B(i-1,u,λ,β,grid))
end
end
```

using

(∂/∂x)f ≈ (1/2) (1/Δx) [F(i+1)-F(i-1)]

(∂/∂x²)f ≈ (1/2) (1/Δx²) [F(i+1)-2⋅F(i)+F(i-1)]

However, my solving doesn’t seem to work at all. For a starter, `sum(sol[i])`

->Inf as `i`

grows. This is plots of the solution for various times:

If I include the final state of the solution it dwarfs everything else:

I solve the PDE using:

```
using DifferentialEquations
u0 = zeros(length(grid)); u0[findfirst(grid.>=10.0)] = 1;
parameters = [10.,0.1]; tend = 10.0;
prob = ODEProblem(cfpe,u0,(0.,tend),parameters);
@time sol = solve(prob,tstops=0:0.1:tend);
```

(The functions are naturally =0 at 0, and I simply assume that they are 0 for very large x (x=200). The solution grows significantly before there’s a significant dnesity at these values, so I don’t think the boundary conditions are a problem)

(The initial condition is everything gathered in a single point. Initially, I thought the discontinuous initial condition might pose a problem, but the solution quickly becomes continuous, and the problem still persists)