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$.
The answer to your question can be a little unintuitive at first so I'll try and explain it as qualitatively as possible.
When you introduce a dopant concentration gradient into a piece of semiconductor, there will initially be a diffusion of carriers from the areas with higher carrier concentration to the areas of lower carrier concentration. Ignoring the diffusion of minority carriers, this results in a charge imbalance as ionised impurities are immobile and cannot diffuse together with the charge carriers.
As a result of this, the ionised impurities attract the charge carriers, establishing an electric field within the semiconductor at thermal equilibrium. Hence the conduction and valence bands are tilted.
However, due to the diffusion of charge carriers, the Fermi level throughout the semiconductor is constant. Due to the tilting of the conduction and valence bands however, its distance from the two bands varies throughout the semiconductor.
Hope this answers your question.
Best Answer
You qualitative reasoning, especially regarding the built-in electric field, is correct but the quasi-fermi level will not simply be upshifted over the whole length of the considered inhomogeneous semiconductor. It will upshifted by $q·V$ at the negative contact. To accurately calculate the current through an inhomogenously doped semiconductor is not a trivial task. In the simplest, one-dimensional, stationary case with a given doping profile $n_0(x)$ , you have to simultaneously solve the nonlinear system of ordinary differential equations for $E$ and $n$ given by:
The drift-diffusion equation $$J=nq\mu_n E+qD_n \frac {dn}{dx}=const \tag 1$$ The Poisson equation $$\frac {d^2 \phi}{dx^2}=-\frac {dE}{dx}= \frac {q (n-n_0)}{\epsilon}$$ Furthermore, you have to assume appropriate boundary conditions at the contacts including the applied bias. In general, such a non-linear system of differential equations can only be solved numerically.