# Numerical solution for a system of PDEs via steady state with time integration method

I'm working to solve the steady-state short circuit current of a solar cells, using the coupled continuity equations with a drift-diffusion expression and Poisson's equation:

D[n[x, t], t] - (1 / q) D[Jn[x, t], x] == G - R; D[p[x, t], t] + (1 / q) D[Jp[x, t], x] == G - R; Jn[x, t] == q [Mu]n (n F + (kb T / q) D[n[x, t], x]); Jp[x, t] == q [Mu]p (p F - (kb T / q) D[p[x, t], x]); D[D[[Phi][x, t], x], x] == (q / ([Epsilon]r [Epsilon]0)) (p[x, t] - n[x, t]);

In particular, I'm following the work from a literature paper that (no surprise) leaves out many details. The paper says to first find the equilibrium solution by solving Poisson's equation numerically with the equilibrium expressions:

nboltz[x] = Nc Exp[(-kb T Log[Nc] eVperJ - q [Phi][x]) / (kb T eVperJ)] ; pboltz [x] = Nv Exp[(-kb T Log[Nv] eVperJ + q [Phi][x] ) / (kb T eVperJ)];

I can get this far using the boundary conditions below (for electric potential). However, it then says: "the steady-state solutions found using the time-evolution method of David's et al, which uses an additional equation for the time derivative of the electric field:

D[F[x, t], t] == - (1 / ([Epsilon]r [Epsilon]0)) (Jn[x, t] + Jp[x, t] - (1 / L) Integrate[ Jn[x, t] + Jp[x, t], {x, 0, L}]);

In addition, I have the boundary conditions,

p[0, t] = Nv n[0, t] = Nc Exp[(Ev - Ec) / (kb T eVperJ)] n[L, t] = Nc p[L, t] = Nv Exp[(Ev - Ec) / (kb T eVperJ)]

Davids et al. says, "To find the steady state solution at an applied voltage bias, a time dependent potential ramp which stops at the desired voltage is applied to the right contact and the equations are integrated forward in time starting from thermal equilibrium until state state is reached . The position independence of the total particle current J = Jn + Jp is sued to verify that steady state has been reached."

Conveniently, I am looking at the short-circuit case where the applied voltage bias is 0. However, I don't understand how one can integrate a steady state solution forward in time. Isn't a steady state solution fundamentally unchanging with time? Moreover, if integrating forward in time is just to help the solver, I am still unsure of how to implement this (I'm working in Mathematica).

Hopefully this all formats well and thanks for the help!

-