DC solution
You can say that you search for the solution in the form of:
$$p(x,t)=p_{dc}(x)+p_{\alpha c}(x,t)$$
So your initial equation will be $$\frac{\partial p_{\alpha c}(x,t) }{\partial t}=D_p\frac{\partial^2(p_{dc}(x)+p_{\alpha c}(x,t))}{\partial^2 x}+\frac{p_n-p_{dc}(x)-p_{\alpha c}(x,t)}{\tau}+\alpha I_0e^{-\alpha x}(1+Me^{i\omega t})$$
then you can decompose into two equations namely they are (you can just add them two see that you have the initial equation):
$$ 0 = D_p\frac{\partial^2p_{dc}}{\partial^2 x}+\frac{p_n-p_{dc}}{\tau}+\alpha I_0e^{-\alpha x}$$ and
$$\frac{\partial p_{\alpha c}(x,t) }{\partial t}=D_p\frac{\partial^2p_{\alpha c}(x,t)}{\partial^2 x}-\frac{p_{\alpha c}(x,t)}{\tau}+\alpha I_0e^{-\alpha x}Me^{i\omega t}$$
Then there are two way of solving the first equation that I see:
Solving the homogeneous equation $ 0 = D_p\frac{\partial^2p_{dc}}{\partial^2 x} - \frac{p_{dc}}{\tau}$ by saying that $p_{dc} = e^{ \lambda x} $, where $\lambda$ gives $L^{-1}$. Next, for a particular solution of a homogeneous equation $0 = D_p\frac{\partial^2p_{dc}}{\partial^2 x}+\frac{p_n-p_{dc}}{\tau}+\alpha I_0e^{-\alpha x}$, you can say that it is $p = p_n - C e^{-\alpha x }$ and get the $C$ constant. So, the solution is the sum of the both homogeneous and non-homogeneous solution.
Other way (without gassing as in 1 case), you can apply Laplace transform over x and you will have less derivatives on x. And the solution will give all the exponents on x
AC solution
As for ac solution I can see the following way:
$$\frac{\partial p_{\alpha c}(x,t) }{\partial t}=D_p\frac{\partial^2p_{\alpha c}(x,t)}{\partial^2 x}-\frac{p_{\alpha c}(x,t)}{\tau}+\alpha I_0e^{-\alpha x}Me^{i\omega t} $$
we apply the separation of variables method $ p_{\alpha c}(x,t) = p_x(x) p_t(t) $, so the equation will be (also I am tired to write a denominator of a partial derivative):
$$p_x \partial_t p_t = D_p p_t \partial^2_x p_x - p_t p_x /\tau + \alpha I_0 e^{-\alpha x } M e^{i \omega t } \, (1)$$
as for $p_x$, I'll say that $p_x = e^{ - \alpha x }$ , and if you apply this solution to the equation (1) and cancelling terms you get
$$ \partial_t p_t = D_p p_t \alpha^2 - p_t /\tau + \alpha I_0 M e^{i \omega t} \,(2) $$
one can use the several model as before, first by
solving homogeneous solution
$$ \partial_t p_t = D_p p_t \alpha^2 - p_t /\tau $$
with guessing $p_t = C e^{\lambda t}$ and getting $\lambda = D_p \alpha^2 - 1/\tau$
solving non-homogeneous equation
we say that simply $p_{ac} = D e^{j \omega t }$ and applying it to the equation (2) we get:
$$ i \omega B = (D_p \alpha^2 - 1/\tau) B + \alpha I_0 M$$
and $B = \frac{\alpha I_0 M }{i \omega - D_p \alpha^2 + 1/\tau}$
So the "ac" solution is
$$ p_{ac}(t) = e^{-\alpha x} ( \frac{\alpha I_0 M }{i \omega - D_p \alpha^2 + 1/\tau} e^{j \omega t } + D e^{(D_p \alpha^2-1/\tau)t }) $$
...
I do not know the dimension of the solutions of differential equation and would appreciate and be curious if someone could correct me if the solution is not exact or there are some mistakes are made.