Diffusion coefficient for asymmetric (biased) random walk
PREFACE
After several edits, this answer provides a naive explanation of why your approach failed, how to fix it (naive-ish) and a completely different (but right) approach to solve the problem.
Intro
You are right: the diffusion coefficient should be $D=4pqD_0$, $D_0$ being the "normal one" (see below for derivations).
I do not know precisely why your approach is not working, but there is a strong evidence something must be wrong which I think depends on the stochastic process involved: if you define ${\Delta x\over \Delta t}=v$ and it is finite, then "your" $D={\Delta x^2 \over 2\Delta t}={\Delta x\over 2} v$ is vanishing in the small $\Delta x$ limit.
The fact is that on the LHS you have a time derivative ($O(\Delta t)$) whereas on the RHS you have a second $x$-derivative ($O(\Delta x^2)$).
Now, since the ratio $v={\Delta x\over \Delta t}$ is finite, it means that $O(\Delta x^2)=O(\Delta t^2)$ so the two sides of the equation do not match. With normal Brownian motion, you would have $\Delta x\sim \sqrt{\Delta t}$ so no problem.
Basically the drift induced by the bias messes things up...
I am disappointed to see that soooo many books and papers do this mistake and use $D_0$ also for biased random walks... (it only makes sense when the bias is external and diffusion is the "normal" one, while in this case it's the same process..)
I now propose you two solutions: the first one is a Fokker-Planck approach. The second one is a Langevin one.
I am sure the results are right (simulations + papers + books prove it) but I am not sure about the steps (as I mainly did them myself). It is hard to find a theoretical treatment of this matter.
Fokker-Planck approach
The problem we are facing here is drift+diffusion: the particle drifts (due to the bias) but also oscillates around its mean value (as normal diffusion would oscillate around $0$).
We define the drift velocity $v={\Delta x\over \Delta t}$. If we had only drift the FP would be: $$\partial_t P_d(x, t)=-v\partial_x P_d$$
Instead if we had only diffusion: $$\partial_t P_D(x, t)=D\partial^2_x P_D$$
Problem: which D do I choose? In general $D\sim \Delta x^2$ but in our case "a fraction" $(p-q)\Delta x$ of the original displacement has been turned to drift, so we renormalize by removing this part (this is equivalent to consider only the variance of the displacement instead that its squared mean value, which again, in normal diffusion are identical so...): $$D={1\over 2\Delta t}\left(\Delta x^2 - (p-q)^2\Delta x^2\right) = 4pq{\Delta x^2\over 2\Delta t}$$ (I used $p+q=1\rightarrow p^2+q^2=1-2pq$).
As you see we found the $D$ I predicted - yet this has been sort of a leap of faith for me, although I see that the diffusion coefficient is sometimes defined as ${var(\Delta x)\over 2\Delta t}={\langle \Delta x^2\rangle - \langle \Delta x\rangle ^2)\over 2\Delta t}$ instead of simply ${\langle\Delta x^2\rangle\over 2\Delta t}$. I never noticed that because in normal diffusion there is no difference as the mean is null.
Anyway, the two processes by themselves are perfectly defined.
We now have to stitch together drift and diffusion. Since the two process happen at the same time, the final probability $P(x, t)$ can be found as:
$$P(x, t)=\int P_D(x-y)P_d(y)dy$$ (meaning: probability of travelling $x-y$ by diffusion times probability of travelling $y$ by drift, integrated over $y$). Again, this is a leap of faith and so is the proof (which I think I managed to find, but wouldn't be sure of it: one just has to take the time derivative of that stuff and the play with the integrals and derivatives a little...) that this leads to:
$$\partial_t P= -v\partial_x P +{D}\partial^2_x P$$ which is the Fokker-Planck we are looking for, with $D=4pqD_0$.
Langevin approach
One of the sources you quoted in a comment already show the following, but I am doing it in a more physical fashion.
I am not sure it can help but I propose you a similar approach (more Langevin-like) which leads to the first two moments of the position $x$ and to an evident definition of the diffusion coefficient.
I am positive about this result.
We suppose a discrete process of $N={t\over \Delta t}$ steps, there $t$ is the total time. At each step, the particles displaces by $\Delta$ with probability $p$ to the right and $q=1-p$ to the left.
Thus at each step: $$x(t+\Delta t)=x(t)+\eta_{\Delta t}$$
where $\eta_{\Delta t}$ is a process which gives $\Delta$ with probability $p$ and $-\Delta$ with probability $q$, such that $\langle \eta_{\Delta t }\rangle = (p-q)\Delta$ and $\eta_{\Delta t}^2=\Delta^2$. Notice that this defines $\eta_{\Delta t}$ only each $\Delta t$ seconds. We do not know what happens at other scales.
So we have, after $\Delta t$: $$dx=x(t+\Delta t)-x(t)=\eta_{\Delta t}$$ and thus its average value is:
$$\langle dx \rangle=\langle \eta_{\Delta t} \rangle=(p-q)\Delta$$ so that $$\langle x(t)-x(0) \rangle=\langle \sum^N_{i=1} x(i\Delta t)-x((i-1)\Delta t)\rangle=\sum_{i=1}^N \langle dx\rangle=N(p-q)\Delta=(p-q){\Delta \over \Delta t}t$$
If we put $v={\Delta \over \Delta t}$ this is your same result and it makes sense. Notice that we had to split $x(t)-x(0)$ in $\Delta t$-sized interval as we do not know what happens on other scales.
Things are harder for the second moment: $$\langle (x(t)-x(0))^2 \rangle=\langle\left(\sum^N_{i=1} x(i\Delta t)-x((i-1)\Delta t)\right)^2\rangle=\langle\left(\sum^N_{i=1} dx\right)^2$$
and that is the square of a sum, so:
$$\langle\left(\sum^N_{i=1} dx\right)^2\rangle=\sum^N_{i=1}\langle dx^2\rangle + 2\sum^N_{i<j}\langle dx\rangle_i\langle dx\rangle_j$$
now, using the moments of $dx$ which we know, being that all moments are equal:
$$\langle (x(t)-x(0))^2 \rangle=N\langle dx^2\rangle+N(N-1)\langle dx \rangle ^2$$ i.e.
$$\langle (x(t)-x(0))^2 \rangle = \Delta^2 N(1+(N-1)(p-q)^2)$$
which can be modified into (using $p+q=1\rightarrow p^2+q^2=1-2pq$):
$$\langle (x(t)-x(0))^2 \rangle = \Delta^2{t\over \Delta t} \left((1-4pq){t\over \Delta t}+4pq)\right) $$
Now, if $p=q=1/2$ you get $$\langle (x(t)-x(0))^2 \rangle = \Delta^2{t\over \Delta t} = 2Dt$$ with $D={\Delta^2\over 2\Delta t}$ which is normal diffusion, as expected.
If instead $p=1$ and $q=0$ (or vice versa)
$$\langle (x(t)-x(0))^2 \rangle = \Delta^2\left({t\over \Delta t}\right)^2 = (vt)^2$$
which again is surely right.
All the intermediate cases are weird. Notice, in addition, that in all other cases if you let $\Delta t, \Delta x\rightarrow 0$ strange things happen.
I guess this is the reson the Fokker-Planck does not turn out right except for simple cases. There must be some "stochastic processes" trick I do now know.
But at least you can rewrite:
$$\langle (x(t)-x(0))^2 \rangle = \Delta^2\left({t\over \Delta t}\right)^2 (1-4pq)+4pq\Delta^2{t\over \Delta t} = (vt)^2(1-4pq)+2D_{eff}t$$ with $D_{eff}=4pq{\Delta^2\over2\Delta t}$.
This process is a drift+diffusion process: the particle drifts with velocity $(p-q)vt$ (which is the position's average value) and around such a value it oscillates. At each time one can compute the variance:
$$\sigma^2=\langle (x(t)-x(0))^2 \rangle -\langle (x(t)-x(0)) \rangle ^2$$ which turns out to be
$$\sigma^2(t) = 2D_{eff}t = 8pqD_0t$$ where $D_0={\Delta^2 \over 2\Delta t}$ is the normal diffusion coefficient.
So I guess the process will be represented by a gaussian distribution (provided that we started with $P(x, 0)=\delta(x-x_0)$ with mean $x0+(p-q)vt$ and variance $\sigma^2(t)$. So the distribution moves and spreads.
Simulations I have been doing right now are in agreement.
Conclusions
I am not sure we found the right Fokker-Planck (or that we found it the right way) but I guess we found $P(x, t)$...! This P(x, t) makes sense and is probably the right one. Yet, as simulations always involve discrete steps, I am not sure about what would happen in the continuous limit though... maybe the validity of some steps may fall apart.
Yet, I think we may consider ourselves satisfied.