# Solving a system of PDE equation in mathematica using a steady-state and time-evolution method [closed]

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 = J_n + J_p$ is used 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!

-

## closed as off topic by dmckee♦Oct 26 '12 at 23:46

Questions on Physics Stack Exchange are expected to relate to physics within the scope defined by the community. Consider editing the question or leaving comments for improvement if you believe the question can be reworded to fit within the scope. Read more about reopening questions here.If this question can be reworded to fit the rules in the help center, please edit the question.

Computational questions are off-topic here, but with your permission I will migrate this to the scientific computation beta site. – dmckee Oct 26 '12 at 23:46
On closer reading, I think this is really two questions, one about what it means to integrate a steady-state solution forward in time, and one about how to implement that integration in Mathematica. The former is, I think, on topic here, and the latter is best placed at Mathematica. So my suggestion would actually be to edit out the Mathematica-specific parts of this question, and then TexanFromTexas, you can post a separate question at Mathematica about the implementation details. What do you think @dmckee? – David Z Oct 27 '12 at 0:24
For what it's worth, I agree that this could be split into a conceptual question about the physics, and a question about the implementation. The former is certainly on topic here, and the latter is likely to do better at Mathematca. – Colin McFaul Oct 27 '12 at 0:29
@DavidZaslavsky I agree. That would be a good way to go. – dmckee Oct 27 '12 at 0:38