As I understand your question, the problem at hand is to solve the PDE for the minority carrier concentration:
$\dfrac{\partial p}{\partial t} = D_p\dfrac{\partial^2 p}{\partial x^2} + \dfrac{p_0 - p}{\tau_p}$
on the half-line $[0, \infty)$, subject to the initial condition:
$p(x, 0) = p_0$
And the boundary conditions:
$p(0, t) = p_0 + N_m$
Where $N_{m}$ is the excess minority carrier concentration due to the photoelectric effect (reference), and
$\lim_{x \to \infty} p(x,t) = p_0.$
To solve this, considering splitting the problem into the homogeneous PDE:
$\dfrac{\partial p}{\partial t} - D_p\dfrac{\partial^2 p}{\partial x^2} + \dfrac{p}{\tau_p} = 0 $
And the inhomogeneous PDE:
$\dfrac{\partial p}{\partial t} - D_p\dfrac{\partial^2 p}{\partial x^2} + \dfrac{p}{\tau_p} = \dfrac{p_0}{\tau_p}$.
Since the system is linear, if we have have a general solution to the homogeneous equation and a particular solution to the inhomogeneous equation, then if their sum satisfies the boundary conditons it will be a solution to the original problem.
It's easy to see that a particular solution to the inhomogeneous equation is $p(x,t)_p = p_0$.
Approaching the homogeneous PDE next, we can use the Laplace transform to find a general solution that satisfies the following BC and IC:
$p(0, t) = N_m$
$p(x, 0) = 0$, and
$\lim_{x \to \infty} p(x,t) = 0$.
Taking the Laplace transform:
$\int_0^{\infty}e^{-st}f(t)dt$
In the time variable of the PDE gives us the ODE:
$s\hat{p}(x) - \mathcal{L}(p(x, 0)) = D_p\dfrac{d^2\hat{p}(x)}{dx^2} - \dfrac{\hat{p}(x)}{\tau_p} \implies $
$\dfrac{d^2\hat{p}(x)}{dx^2} = \hat{p}(x)\left(\dfrac{s}{D_p} + \dfrac{1}{D_p \tau_p}\right) = $
$\dfrac{d^2\hat{p}(x)}{dx^2} = C_n\hat{p}(x)$, where
$ C_n = \left(\dfrac{s}{D_p} + \dfrac{1}{D_p \tau_p}\right)$.
The general solution to this equation can be found by standard ODE methods and has the form:
$\hat{p}(x) = C_1e^{\sqrt{C_n}x} + C_2e^{-\sqrt{C_n}x}$.
Since we only desire solutions that decay as $x \to \infty$, the first solution is nonphysical and $C_1 = 0$.
The initial condition $p(0) = N_m\implies \hat{p}(0) = \dfrac{N_m}{s} \implies C_2 = \dfrac{N_m}{s}$,
so the solution to the ODE in the S domain is:
$\hat{p}(x) = \dfrac{N_m}{s}e^{-\sqrt{C_n}x}.$
Approaching the right hand side first, with some help from Wolfram Alpha, we find that the inverse Laplace transform of the function $P(s) = \frac{A}{s}e^{-\sqrt{s + \alpha}{x}}$ is:
$\frac{1}{2}Ae^{-\sqrt{\alpha}x}\left(\mathrm{erf}\left(\frac{2\sqrt{\alpha}t - x}{2\sqrt{t}}\right)+ e^{2\sqrt{\alpha}x}\mathrm{erfc}\left(\frac{2\sqrt{\alpha}t + x}{2\sqrt{t}}\right) + 1\right)$ for $x > 0$,
where $\mathrm{erfc}(z) = 1 - \mathrm{erf}(z)$ is the complementary error function.
Our function $\hat{p}(x)$ is of the form $\frac{A}{\beta}\frac{1}{{\frac{s}{\beta}}}e^{-\sqrt{\frac{s}{\beta} + \alpha}{x}}$, so applying the scaling theorem:
$\mathcal{L^{-1}}(P(\frac{s}{\beta})) = \beta p(\beta t)$ we finally get:
$\mathcal{L^{-1}}(\hat{p}(x)) = p(x,t)_h = \frac{1}{2}N_me^{-\sqrt{\alpha}x}\left(\mathrm{erf}\left(\frac{2\sqrt{\alpha}t\beta - x}{2\sqrt{t\beta}}\right)+ e^{2\sqrt{\alpha}x}\mathrm{erfc}\left(\frac{2\sqrt{\alpha}t\beta + x}{2\sqrt{t\beta}}\right) + 1\right)$,
where $\alpha = \frac{1}{D_p\tau_p}$, $\beta = D_p$, and the subscript "h" in $p(x,t)_h$ indicates that this is the homogeneous solution.
Summing the particular and homogeneous solution we get:
$p(x,t)_h + p(x,t)_p = p(x,t) = \frac{1}{2}N_me^{-\sqrt{\alpha}x}\left(\mathrm{erf}\left(\frac{2\sqrt{\alpha}t\beta - x}{2\sqrt{t\beta}}\right)+ e^{2\sqrt{\alpha}x}\mathrm{erfc}\left(\frac{2\sqrt{\alpha}t\beta + x}{2\sqrt{t\beta}}\right) + 1\right) + p_0$,
which satisfies the conditions of the original problem. It should be straightforward to take the limit of this expression as $t \to \infty$.
Best Answer
The answer to this paradox is, as far as I can see, the following. The so called ambipolar transport equation is derived in Neaman's book by assuming (as an approximation) strict charge neutrality so that the local excess carrier densities for electrons and holes become exactly equal. Let us consider the relevant case of a p-type semiconductor with low electron injection (minority carriers). If you have a local excess electron density the assumption means that you have exactly the same excess hole density in order to keep neutrality. As the ambipolar diffusion is dominated by the excess electron density giving the diffusion coefficient for electrons in this approximation you get exactly the same spatial and temporal behavior of the excess hole concentration. You cannot use the diffusion current density equation for holes anymore because it contradicts the strict neutrality assumption. In reality you have only quasi-neutrality, the spatial shape of excess carrier densities is not exactly the same so that a space charge is created leading to internal fields so that the hole (and electron) current has also a drift component in addition to the diffusion component.