The system of coupled equations describing process of electrodiffusion (\ref{electrodiffusion}) written in terms of the FEM method (see Eq.~(\ref{pnp_final})) where k_{+} = k_{-} = k, D_{+} = D_{-} = D and the time step equals 0.01 has been numerically solved up to the final time T_{end} by using the Newton algorithm. The parameter k denotes (ke)/(ε_{0}ε) where k = (De)/(k_{B}T). Computations have been performed for below-given sets of boundary conditions and parameters. The first considered case shows boundary layers in x, y - directions. The second one presents the range of k values in which the one stable solution to this nonlinear problem has been obtained for the assumed values of D = 0.05 and Δt = 0.01. Moreover, solution's convergence depends on h_{0} - the mesh parameter and the order of FEM approximation. The third case is computed for the physical values of parameters which locate solutions within the stable regime.
In all cases an initial guess of n_{+}, n_{-} and φ distributions has been chosen as 0 everywhere in the domain apart from its boundaries.
Fig. 7 presents obtained profiles of cations n_{+} and the potential φ at the center of the domain i. e. at z = π/2 and at the final time T_{end} computed for the first and the third case. There is no visible difference between cations n_{+} and anions n_{-} distributions in both cases so the latter is not shown.
In the first case, the maximum of n_{+, i+1} - n_{+, i} is decaying from 0.023 for i = 0 (t = 0) to 0.0093 for i = 38 (t = 0.38) and in the latter from 0.17 (✗ 10^{7} [n.u.]) at t = 0 to 0.0033 (✗ 10^{7} [n.u.]) for i = 398 (t = 3.98 s). The system of particles tends to its stationary state. Moreover, the maximum of difference between n_{+} and n_{-} distributions computed at each node at the time T_{end} equals 1.3e-09 in the first case and 2.3e-10 (✗ 10^{7} [n.u.]) in the latter one. Employing physical notion it means that the system of charged particles is electroneutral.
The picture shows FEM approximation (linear) of the system of electrodiffusion equations being the boundary value problem. The solutions of a) and c) the distribution of cations n_{+} and b) and d) the profile of the potential φ are depicted in the center of the cubic domain i. e. at z = π/2 and at times T_{end} = 0.39 with the values of parameters Δt = 0.01 [units], k_{+} = k_{-} = 0.05 [units], D_{+} = D_{-} = 0.05 [units] (upper cases) and the physical ones (lower cases), respectively; computations have been performed on the uniform mesh with the volume of tetrahedron V_{0} = 0.01 and the element size h_{0} = 0.44 (upper cases) and V_{0} = 0.005 for h_{0} = 0.35 (lower ones).