A fast and robust numerical scheme for solving models of charge carrier transport and ion vacancy motion in perovskite solar cells

Drift-diffusion models that account for the motion of both electronic and ionic charges are important tools for explaining the hysteretic behaviour and guiding the development of metal halide perovskite solar cells. Furnishing numerical solutions to such models for realistic operating conditions is challenging owing to the extreme values of some of the parameters. In particular, those characterising (i) the short Debye lengths (giving rise to rapid changes in the solutions across narrow layers), (ii) the relatively large potential differences across devices and (iii) the disparity in timescales between the motion of the electronic and ionic species give rise to significant stiffness. We present a finite difference scheme with an adaptive time step that is posed on a non-uniform staggered grid that provides second order accuracy in the mesh spacing. The method is able to cope with the stiffness of the system for realistic parameters values whilst providing high accuracy and maintaining modest computational costs. For example, a transient sweep of a current-voltage curve can be computed in only a few minutes on a standard desktop computer.


Introduction
Recent rapid improvements in power conversion efficiency have brought metal halide perovskite solar cells (PSCs) to the forefront of the emerging thin-film photovoltaic technologies. Efficiencies in excess of 20% [4,32], comparable to that of standard crystalline silicon devices, have been achieved using methylammonium lead tri-iodide and other perovskite materials as absorber layers in thin film architectures [15,17]. This high performance, together with their relatively low cost of manufacture, mean that the continued development of PSCs is an extremely active area of research. Recent reviews of the field have been given in [20,22,24,34,35]. Planar PSC architectures are formed by a perovskite absorber layer sandwiched between an electron transporting layer (ETL) and a hole transport layer (HTL). The most commonly used perovskite absorbing material is methylammonium lead tri-iodide (CH 3 NH 3 PbI 3 ) [14], however, recently other mixed halide formulations (in which some of the iodide ions are substituted by other halides) and mixed cation formulations (in which the methylammonium cation is fully or partially replaced with formamidinium and/or caesium) have also been successfully employed [25,29,32]. Numerous different materials have been used as ETLs and HTLs but common choices are titania (TiO 2 ) for the former and spiro (2,2',7,7'-tetrakis-(N,N-di-p-methoxyphenylamine)-9,9'-spirobifluorene) for the latter.
Under illumination, incident photons with energies above the band gap are absorbed in the perovskite absorbing layer generating weakly bound excitons (binding energy ∼50meV [16]). These excitons readily dissociate into a free electron in the conduction band and a hole in the valence band which can move independently under the influences of both thermally-induced diffusion and electronically-induced drift. Since (i) the conduction band in the HTL is above that in the perovskite and (ii) the valence band in the ETL is below that in the perovskite, barriers exist for the entry of holes and electrons into the ETL and HTL respectively. The differences in the band energies of the different semiconductors also lead to the formation of a built-in field in the perovskite which encourages electrons to move towards the ETL and holes towards the HTL, thereby setting up a useful electronic current. A cartoon of this process is shown in figure 1.
One notable peculiarity of PSCs is their long timescale transient behaviour occurring on the order of 10s of seconds. This is observed in (i) current-voltage (J-V) curve 'hysteresis' [31,36] and (ii) current transients [23,36]. Initially, three possible explanations for these slow dynamics were postulated [31], namely: (a) the formation of ferroelectric domains (b) large-scale trapping of electrons, and (c) mobile ions. However, as discussed in [27], there is a growing consensus that the only mechanism that offers a coherent explanation of the observed data is the slow motion of positively charged anion vacancies. More recently, very long timescale reversible transients in the efficiency of PSCs have been observed over periods of several hours [6]. These have been attributed to the motion of slow cation vacancies.
A range of approaches exist to modelling PSC behaviour that extend from fundamental atomistic density functional theory (DFT) calculations (e.g. [8]) to equivalent circuit device models (e.g. [3]). High computational costs mean that DFT can only be applied to a few atoms and over extremely short timescales. At the other extreme, equivalent circuit models, which are very straightforward to solve, are hard to connect directly to the device physics. An intermediate path, which leads to a computationally tractable model that can be directly interpreted in terms of device physics, is given by charge transport, or equivalently driftdiffusion, modelling. Notably this approach allows (as in [27]) for parameters obtained from DFT calculations on the perovskite structure to be incorporated into the model. It has also been used in a variety of other solar cell applications (for example in organics [1,9,28]). Initially the importance of ion motion in PSCs was not fully appreciated and as a consequence charge transport models, that treated only electron and hole transport, were formulated [10,18]. Subsequently models that considered the coupling between ion vacancy motion and charge transport were investigated in [2,6,12,13,23,27,37]. However the incorporation of realistic (high) densities of ion vacancies leads to a model that is computationally challenging to solve owing to (i) narrow Debye layers across which rapid changes in solution occur, (ii) very large changes in the magnitude of the solution across the device and (iii) the large disparity between the timescales for ion vacancy motion and electron and hole transport. Of the aforementioned works on charge transport modelling of PSCs, only [5,23,27] manage to obtain solutions in physically relevant parameter regimes. However in all these studies the analysis relies on making some level of asymptotic approximation to the governing equations. Whilst these approximations are well-justified in most scenarios, and lead to accurate approximate solutions, a numerical treatment of the full system of equations, that is capable of furnishing solutions across all relevant timescales and operating regimes without simplification, is highly desirable. This motivates the aim of this work which is to present a numerical method capable of accurately solving such problems using modest computational power.
The work is set out as follows. In the next section, we outline the charge transport model in dimensionless form, describe how the system is related to its dimensional counterpart and give estimates for the sizes of the dimensionless parameters. In §3 we describe the spatial (finite difference) discretisation. Then, in §4, we outline implementation of the method in MATLAB [19] using the Advanpix Multi-Precision Toolbox [21]. In §5 we benchmark the scheme using some experimentally motivated test cases, comparing the accuracy of the numerical approach to the asymptotic approach developed in [5,26] and demonstrating the pointwise convergence of the method. In §6, we briefly discuss how this approach compares to other methods that we tried, including spectral methods and MATLAB's built-in solver PDEPE. Finally, in §7, we present our conclusions.

The model
Here we follow [27] and consider a model that has been formulated for charge transport in a planar PSC consisting of a perovskite absorber layer sandwiched between highly doped electron and hole transport layers. The doping, and the resulting high conductivity of the transport layers, allows us to treat these layers as 'quasi-metals' through which the electric potential is uniform and equal to that on their respective contacts. In turn this leads to a single-layer model in which all the relevant physics takes place within the perovskite layer. The key physical processes described by the model are the motion, generation and recombination of highly mobile charge carriers (electrons and holes) and their interaction with (much less mobile) anion vacancies and a uniformly distributed stationary cation vacancy distribution.
Since the focus of this work is on the numerical method required to solve this model, and detailed descriptions are given in [27,5], the model is only outlined here. Notably, the numerical scheme that we develop for the single-layer model can easily be extended to more realistic multi-layer models. For example, a three-layer that incorporates charge carrier dynamics in the electron and hole transport layers. An investigation of this three-layer model will be the subject of future work.

Model equations
In order to present the method in the simplest fashion possible, whilst highlighting potential numerical difficulties, the model equations are presented in dimensionless form (for full details of the non-dimensionalisation see [5]). The key variables in the problem are time, t; the perpendicular distance from the interface of the perovskite layer with the ETL, x; the mobile anion vacancy density, P ; the free-electron density, n; the hole density, p; the electric field (in the x-direction), E; the electric potential, φ; the anion vacancy flux (in the x-direction), F P ; and the electron and hole current densities, j n and j p , respectively. The (dimensionless) model consists of three conservation equations for the three mobile charged species: in which G(x) and R(n, p) represent rates of charge carrier generation and recombination, respectively. These conservation equations couple to Poisson's equation: The term (P − 1 + δ (p − n)), in the equation above, is the dimensionless charge density and comprises contributions from anion vacancies, stationary cation vacancies, holes and electrons, respectively. Equations (1)-(4) hold within the perovskite layer which is bounded by its interface with the ETL, on x = 0, and its interface with the HTL, on x = 1. In practice, the three dimensionless parameters λ (the ratio of the Debye length L d to the the perovskite layer width), ν (the ratio of the timescales for electron and hole motion to that for ion vacancy motion) and δ (the ratio of typical carrier concentration to typical vacancy concentration) are all very small, see (11).

Boundary and initial conditions
At the interfaces with the ETL (on x = 0) and the HTL (on x = 1) we require that there is no flux of anion vacancies and that the potential is specified (and equal to that in the adjacent contact). In additionn, the free electron concentration on x = 0, is determined by the band energy offsets with the ETL; whilep, the hole concentration on x = 1, is determined by the band energy offsets with the HTL; and j p | x=0 and j n | x=1 are determined by R l and R r , the rates of surface recombination on these interfaces. It follows that the boundary conditions read: Here the total potential drop across the cell, Φ bi − Φ(t), is split into two parts; a built-in voltage, Φ bi , that arises from the difference in band energies between the ETL and the HTL and an applied voltage, Φ(t). In this formulation of the problem, potentials are measured in units of the thermal voltage (approximately 0.025V). The built-in voltage is typically around 1.1V and standard experiments vary the applied voltage within the range of −0.5 to 2V, therefore the total potential drop across the cell, Φ bi − Φ(t), can be quite large. In turn this can lead to very large variations in both ion vacancy and carrier concentrations across the perovskite layer which renders furnishing numerical solutions challenging. The problem is closed by the following initial conditions for the electron, hole and anion vacancy concentrations:

Definition of dimensionless variables and parameters
The dimensional counterparts of the variables used in the model (1)-(7), which we denote by a star, are retrieved from the following rescaling (see [5]): Here b denotes the width of the perovskite layer, L d the Debye length, D + the anion vacancy diffusion coefficient, q the elementary charge, V T the thermal voltage (i.e. k B T /q where T the absolute temperature), F ph the incident photon flux,D a typical electronic charge carrier diffusivity and N 0 the equilibrium anion vacancy density. The Debye length is calculated based on the most populous charged species, which in this instance is the anion vacancies, as follows: where ε is the permittivity of the perovskite. The dimensionless parameters in (1)- (7) are given in terms of physical constants by the relations: Here D n and D p are free-electron and hole diffusivities, respectively, and n b and p b are the (dimensional) electron density on the ETL interface and the (dimensional) hole density on the HTL interface, respectively. In [5] it is shown that, for a typical planar device formed by a methylammonium lead tri-iodide perovskite absorber layer sandwiched between a titania ETL layer and a spiro HTL layer, the dimensionless parameters defined in (10) take the values The difficulty in solving (1)-(3) arises from the extreme values of the parameters ν, λ and Φ bi − Φ(t). The very small value of ν reflects the large disparity in timescales for the electronic (fast) and ionic (slow) motion. This feature of the problem necessitates the use of an adaptive timestep since any fixed time stepping method would either be prohibitively slow or be incapable of capturing the fast electronic dynamics. As is typical for many electrochemical problems (see e.g. [11]), the parameter characterising the relative size of the Debye length, in this case λ, is also extremely small. This gives rise to appreciable stiffness in the system owing to the rapid changes in the solution across the narrow Debye layers. This issue is further exacerbated by the relatively large value of Φ bi − Φ(t). The concentrations of the different charge species in the slender Debye layers are approximately Boltzmann distributed, i.e. they vary exponentially with the potential. Thus, a change in the potential of size 20 (the size of the drop expected across one Debye layer when the drop across the whole device is 40, e.g., near short-circuit) gives rise to a change in n, p or P by a factor of exp(20) = O(10 9 ) across a region of width O(10 −3 ). Therefore, what at first glance appears to be merely a stiff problem is actually an extremely stiff problem. Notably,these very large changes in the magnitudes of the concentrations of the charged species can give rise to large condition numbers and, depending on the computational precision being used, significant round off errors.

Numerical scheme
The central technique underlying our numerical scheme is the method of lines. The spatial derivatives in (1)-(7) are discretised using second order accurate finite differences on a 'staggered grid'. Since the equation governing the electric potential is elliptic, in contrast to those for the charge carrier densities which are parabolic, the discrete system takes the form of a system of coupled differential-algebraic equations (DAEs). As such, it requires a specialised algorithm for temporal integration. Here, we employ MATLAB's integrator ode15s [19] to evolve in time. In the physically relevant parameter regimes, typical solutions exhibit rapid changes in the narrow Debye (boundary) layers which gives rise to significant stiffness in the DAE system. This difficulty is overcome by (i) employing a non-uniform grid spacing, and (ii) using Advanpix's Multiprecision Computing Toolbox [21] to overload MATLAB's native commands with arbitrary-precision counterparts, thereby retaining good accuracy despite the large condition numbers of the underlying matrices. We note that the treatment of the conservation equations used here is in contrast to the Scharfetter-Gummel scheme that is widely used in electrochemical problems [30]. This difference arises because their discretisation is aimed at addressing issues of transport within drift-diffusion equations whereas our is dealing with the difficulties associated with stiffness.

Computational grid
The computational grid is comprised of N + 1 arbitrarily positioned grid points, denoted by x = x i for i = 0, ..., N, which partition the domain x ∈ [0, 1] into N subintervals. We further introduce a set of N 'half-points' denoted by x = x i+1/2 and defined as follows A sketch of the grid is shown in figure 2. The anion vacancy profiles are determined subject to the Neumann conditions (5a) and (6a) which require zero flux of ions at each boundary. Thus, it is natural to compute the ion vacancy flux at the grid points and the ion vacancy density at the half points. Not only does this allow (5a) and (6a) to be imposed straightforwardly, it also ensures that the property of global conservation of anion vacancies is inherited by the discrete system 1 . The equation determining the electric potential is to be solved subject to the Dirichlet conditions (5b) and (6b) which motivates tracking the potential at the grid points and the electric field at the half points. Since the governing equations for electrons (holes) are to be solved subject to one Dirichlet and one Neumann condition, namely (5c) and (6d) ((5d) and (6c)), it is not so clear whether the concentration or flux should be defined at the grid points. Here, we elect to track the concentration at the grid points and the fluxes on the half-points so that evaluating the electron and hole contributions to the charge density on the right-hand side of (4a) can be done so without interpolation, thereby  avoiding additional errors. In summary we introduce discretised variables defined by

Spatial discretisation
Spatial discretisation will be carried out by introducing the discrete operators D i and D i+1/2 which approximate the first derivative operator at x i and x i+1/2 , respectively, and the discrete averaging/interpolation operators, I i and I i+1/2 which approximate half points quantities on the grid points and grid point quantities at the half-points, respectively. For a generic variable w, these four operators are defined as follows: A standard error analysis indicates that each operator is second order, i.e., that the errors are expected to decrease proportional to the grid spacing squared. The discrete operators defined by (17)- (20) can be used to approximate the electric field, and electron and hole currents at the half-points. The discretised versions of equations (4b), (2b) and (3b) are The anion vacancy flux can then be computed on both the internal and boundary grid points by discretising equation (1b) and its boundary conditions (5a) and (6a). We have The ODEs governing the evolution of the anion vacancy density, arising from (1a), are for i = 0, ..., N − 1.
The algebraic equations for the potential resulting from discretisation of (4a) and the boundary conditions (5b) and (6b) are It is straightforward to discretise (2a) and (3a) to formulate ODEs governing the evolution of the electron and hole densities on the internal grid points. For the electrons (holes) a Dirichlet-type boundary condition is to be imposed on x = 0 (x = 1), see (5c) ((6c)), whilst a specified value of the flux is to be imposed at the boundary x = 1 (x = 0), see (6d) ((5d)). The former of these conditions immediately specifies n 0 =n (p N =p) whilst the latter is imposed by linear extrapolation of the electron (hole) current to the boundary x = 0 (x = 1). Hence we arrive at the following system of ODEs and algebraic equations which are sufficient to determine the values of electron density at all grid points.
for i = 1, .., N − 1, Similarly, the equations determining the hole density at each grid point are for i = 1, .., N − 1, We note that the boundary conditions (32) and (33) are the only conditions not imposed exactly. However, the extrapolation involved only introduces second order errors, and therefore the scheme as a whole is still second order as demonstrated by figure 4.

Implementation
The system of DAEs formulated above is evolved forward in time using MATLAB's ode15s.
In order to minimise computational cost, by minimising the size of the system, we eliminate superfluous variables (namely E i+1/2 , j p i+1/2 , j n i+1/2 and F P i ) and assemble the remaining 4N + 3 unknown functions of time into one column vector u(t) as follows: where a superscript T denotes a transpose. In (37) P is a column vector of length N whose entries are the functions of time P i+1/2 (t) for i = 0, ..., N −1. Similarly, Φ, n and p are column vectors of length N + 1 whose entries are the functions φ i (t), n i (t) and p i (t) respectively. The problem to be solved can now be written in the form where f(u) is a nonlinear vector function of length 4N + 3 whose first N entries are the right-hand sides of (26), the subsequent N + 1 entries are the right-hand sides of (27)-(29), followed by the N + 1 right-hand sides of (30)- (32), and finally the N + 1 right-hand sides of (33)- (35). The matrix M is a (4N + 3) × (4N + 3) diagonal mass matrix whose entries are the coefficients of the time derivative terms in the equations (26)-(35), see figure 3. Since the governing equations for the values of φ i , (27)- (29), and the discrete boundary conditions (30), (32), (33) and (35) contain no temporal derivatives the corresponding entries on the diagonal of M are zero and hence the mass matrix is singular. It is this feature of the system that renders the problem a system of DAEs and motivates our choice of solver, namely ode15s, which is one of relatively few solvers that can handle problems of this type. Moreover, it offers an adaptive timestep which is able to deal with the disparity in timescales between the electronic and ionic motion both with good accuracy and without high computational penalties. The ode15s time step requires numerical approximation of the Jacobian of f. However, it is clear from the structure of the discrete system that many entries in the Jacobian are zero. A heavy reduction in computational effort, around 0.005N 3/2 in the convergence test case, can be achieved by exploiting the jpattern option which allows the user to 'flag' only a subset of the entries in the Jacobian matrix that need to be numerically approximated (those that are not flagged are assumed to be zero) at each time step, see figure 3.
When simulating many of the regimes of physical interest, we find that the condition number of the Jacobian is large, sometimes O(10 16 ) or greater, which has the potential to severely hamper the accuracy of matrix inversions performed by the solver, ode15s. To overcome this difficulty we make use of a third-party toolbox, Advanpix, which extends MATLAB's functionality to work with floating point numbers of higher precision. This can result in a loss in the speed of computation for simple cases (those which can be computed accurately on a small number of grid points) for which runtimes, originally of a few seconds on a standard desktop computer, are typically 15-30 times longer. However, the higher precision offered by Advanpix's Multiprecision Computing Toolbox [21] is crucial in allowing the numerical scheme to cope with the stiff regimes of practical interest.

Verification and test cases
We verify the performance of the scheme using physically relevant test cases. Here we choose to compare our numerical results to the asymptotic solution given in [5] in the particular case where the bulk recombination, R(n, p), is monomolecular (depending solely on the local hole concentration p) and the surface recombination rates, R l (p) and R r (n) are both zero, such that Here γ is a dimensionless rate constant. As noted in [5], monomolecular hole-dominated bulk recombination is the limit of the Shockley-Read-Hall (SRH) recombination law when electron pseudo-lifetime is much less than hole pseudo-lieftime. This is a reasonable (though not completely accurate) description of recombination in the perovskite material methylammonium lead tri-iodide [33]. Note however that the numerical scheme presented here is readily able to deal with nonlinear recombination rates such as the full SRH law. As is typical, in such applications, we assume that the photo-generation rate, G(x), obeys the Beer-Lambert law, which has the dimensionless form in which 1/Υ is the dimensionless absorption length. Estimates of these from [5] are A simple choice of initial conditions satisfying the boundary conditions is In practice, the choice of initial conditions is not particularly important for the modelling of PSCs as, in any experimental procedure, it is standard practice to include a pre-conditioning step. This step involves holding the applied potential Φ(t) constant for a sufficiently long time such that any initial transients associated with charge carrier and ion vacancy motion decay towards zero.

Choice of spatial grid
In anticipation of the fact that the largest gradients in the solution appear in narrow Debye layers adjacent to the domain boundaries, we use a spatial grid comprised of Chebyshev nodes on the interval [0, 1] defined as follows By concentrating points near the domain boundaries we are able to achieve good resolution in the Debye layers without wasting computational effort by over resolving in the bulk where the solution varies more slowly.

Spatial convergence
Here we verify the expected second order pointwise convergence of the scheme by testing the method in a scenario in which an illuminated cell is initially operating at an applied voltage equal to the built-in voltage Φ = Φ bi = 40 (see 11), so that the potential drop across the perovskite layer is zero. The applied voltage is then rapidly decreased to Φ = 20 so that the device is running in its power generating regime. We accomplish this by defining Here β characterises the timescale for altering the applied voltage and t end is the time at which the simulation is terminated. We terminate the simulation at t = t end which is not so large that the solution has reached steady state. During this convergence test we are therefore examining the scheme's accuracy during the transient and not just at steady state.
These are chosen because, independent of the number of grid points used, they are available immediately following the temporal integration without the need for an additional interpolation step, thereby allowing a direct interrogation of the scheme. Due to the lack of exact solutions to the model, it is necessary to verify convergence by measuring the error relative to a solution computed on a highly refined grid. For some scalar quantity v (which could be any of those defined in (45)) computed using the scheme on N grid points, denoted by v (N ) , we can define the error, for some N M ≫ N. The pointwise convergence of the quantities in (45) is shown in figure 4. These results which demonstrate second order order convergence of the scheme, as predicted by the naive analysis of (17)- (20).

A current decay transient: comparison between numerical and asymptotic solutions
Here we simulate a cell that is preconditioned at Φ = Φ bi for a sufficient time to eliminate transients, before undergoing a smooth but rapid decrease in applied bias from Φ = Φ bi at t = 0 to Φ = 0 some short time later, obeying, with β = 10 6 , t end = 1.
A comparison between the photocurrent calculated from the numerical solution (as described above) and the asymptotic solution (as described in [5]) is made in figure 5, showing remarkable agreement between the methods for all but very small times. The inset shows the very short time behaviour where there is significant discrepancy between the two solutions. This results from the numerical method capturing fast electronic transients (associated with electron and hole motion) that are absent from the asymptotic solution. The ability to simulate on the timescale of electronic transport through to that of ion motion is important for rigorous investigation of many properties of PSCs. In figure 6, we demonstrate good agreement between the numerical and asymptotic solutions, with the noticeable exception of panel (c) for the electron distribution in which the solutions vary in magnitude, but not shape. From a convergence check (increasing both spatial and temporal accuracy), it can be confirmed that the numerical scheme has indeed converged. The discrepancy arises from a small O(10 −2 ) error in the asymptotic solution to the electric potential which in turn leads to the significant error in the asymptotic solution  for electron density n, that can be seen panel (c).

A current-voltage curve: comparison between numerical and asymptotic solutions
Measurements of the current generated by a solar cell in response to a backwards and forwards sweep of the applied voltage are a common way of measuring cell properties and, in the case of PSCs, frequently result in significant sweep-rate dependent 'hysteresis' [31]. This 'hysteresis' is a signature of the slow timescale ion motion in the cell. The standard experimental procedure for generating such a J-V hysteresis-curve, including a preconditioning step, will be simulated as follows. The cell is preconditioned in forward bias (at 1.2V, corresponding to Φ ≈ 46.7), then the voltage is scanned smoothly from forward bias (Φ ap > Φ bi ) to short circuit (Φ ap = 0) and back at a scan rate equivalent to 100 mV/s (corresponding to a rate of ≈ 7.1 in dimensionless units). The resulting J-V curve is compared to the corresponding asymptotic solution from [5] in figure 7. Plots of the solutions variables in space, at equally-spaced times through the hysteresis sweep, are made in figure 8.

Discussion
We previously attempted to solve this stiff numerical problem using spectral methods via the Chebfun package [7] for MATLAB. An advantage of spectral methods is the high spatial accuracy that can be obtained with a small number of grid points. However, this method requires such fine resolution in time to maintain accuracy that it becomes too computationally expensive to solve problems in the desired parameter regimes. In particular, the ratio of the Debye length to the perovskite layer width, λ, caused a significant slowdown if taken to be O(10 −2 ) or smaller, and as a result solutions for the realistic value of λ = 2.4 × 10 −3 were unobtainable in a reasonable amount of time (days). Calado et al. [2] use MATLAB's built-in solver PDEPE to investigate a drift-diffusion model for PSCs. However, we found that in realistic parameter regimes this method has a condition number that is so large that matrix inversions cannot be performed to within even single-digit accuracy. The method could thus only be applied to toy problems in order to give qualitative indications of the physical phenomena occurring in real cells. The increases in accuracy afforded by the staggered grid and the increased arithmetic precision used in our approach are essential in order to investigate the behaviour of real PSCs. A secondary disadvantage of using PDEPE is that it can only accept certain types of boundary conditions and this prevents the extension of the method to a three-layer model that incorporates charge carrier transport in the ETL and the HTL of a PSC. In particular PDEPE is unable to deal with conditions posed on internal boundaries such as occur at the ETL/perovskite and the perovskite/HTL interfaces. In contrast, the tailored numerical scheme detailed in this work can be readily extended to include more layers and incorporate any necessary type of interface condition.

Conclusions
We have presented a numerical scheme for solving a PDE drift-diffusion model of ion vacancy and charge carrier transport in a perovskite solar cell. The scheme uses second-order finite difference approximations of spatial derivatives, on a non-uniform staggered mesh, to reduce the system of PDEs to a (large) system of ordinary differential-algebraic equations (with time as the independent variable). This system of ordinary differential-algebraic equations is solved with the MATLAB routine ode15s in conjunction with Advanpix's Multiprecision Computing Toolbox [19,21].
For realistic parameter values and in appropriate operating regimes, the problem exhibits significant stiffness owing to (i) the small Debye length, (ii) large potential differences across the device (giving rise to large and rapid changes in solution across narrow Debye layers), and (iii) vastly different timescales for the transport of ion vacancies and electronic charge carriers. We were able to circumvent these difficulties by (a) using a non-uniform mesh to selectively concentrate grids points in the Debye layers in order to obtain high accuracy without prohibitive computational cost, and (b) using an adaptive timestep. In doing so we have presented the only published method for numerical solution of a truly realistic charge transport model for a metal halide perovskite solar cell. Moreover, the method can be used to simulate experimental protocols, such as current decay transients and current-voltage sweeps, in only a few minutes on a desktop computer.