BVP with sharply varying values

Dear Matlab users,

I am currently trying to solve the drift-diffusion equations in a n-i-p structure. The domain is 1D and divided in 3 parts:

• n region from x=0 to xn
• i region from x=xn to xper
• p region from x=xper to 1

A simplified version of the equations on the functions n(x), p(x) and V(x) can be given as such (scaled with as few parameters as possible to reproduce the problem):

I am solving them as a BVP problem using bvp5c solver (more on boundary conditions below).

Important remarks on the parameters:

• lambda^2 has a small value (1e-5)
• G and Ndop are piecewise constant. They vary very sharply (discontinuously) at the interfaces of the n-i-p structure. And i think the problem boils down to this.

Here is my my function implementing the equation:

function dydx = diffeq(x,y)
 % the function vector is as follows : y = [p; dp/dx; n; dn/dx; V; dV/dx] 
 if x<xn %if we are in the n region
 Ndop=Na;
 G=0;
 elseif x<xper %if we are in the i region
 Ndop=0;
 G =Gcons;
 else %if we are in the p region
 Ndop=-Na;
 G=0;
 end
 d2V=(y(3)-y(1)-Ndop)/lamb2;
 dydx=ones(length(y),1);
 dydx(1) = y(2);
 dydx(2) = ( -G - y(2)*y(6) - y(1)*d2V ) ; 
 dydx(3) = y(4);
 dydx(4) = ( -G + y(4)*y(6) + y(3)*d2V ) ;
 dydx(5) = y(6);
 dydx(6) = d2V;
end

G is 0 in n and p regions, and 8e-5 in i region. This discontinuity used to cause problems, but since I re-scaled the equations to have more resonnable numerical values, this issue seemsfixed. Ndop is a constant Na in the n region, 0 in the i region and -Na in the p region. This works ok for low values of Na (up to ~1e-3), but higher values cause the convergence to fail with messages like "singular jacobian encountered" or "Unable to meet the tolerance without using more than 2000 mesh points" (depending on the exact value). I have tried increasing maximum mesh size, or to relax tolerance, but it does not solve the issue.

I would eventually need a Na value of 1 so that Ndop varies as shown on the picture below

Question: Am I doing something wrong in the use of bvp5c? Is there a way to fix this issue so that the solver can handle the discontinuity, or will it be necessary to implement more specialized solving algorithm for those equations?

The boundary values that I am using are given below. (NB: they are not the most common assumptions for this system, but I have seen them used in a paper and wanted to give it a shot. The most "canonical" BC do not work any better anyhow):

I have attached my matlab files for testing. The file "main" will call the "solver" file 2 times; one with low Na to show convergence, and one with higher values causing the problem. the second run uses solution of the first one as initial guess. A temporary file is used to contain the parameters. I have implemented the jacobians (I am not 100% sure that they are fully accurate, I will check further), but the problem shows up regardless of whether I use them or not so the root cause is something else.

[David Wilson]

If the problems really are due to the sharp disontinuity in Na, then one approach would be to soften it with say a tanh(kx) function that smooths your step discontinuity. You get to adjust the steepness with the parameter k. Start with k small, and gradually increase it. Use the previousl solution from bvp4c (perhaps not 5c) as the initial guess for the new solution. Sneak up on the solution.