How to solve continuity equations together with Poisson equation?

As I'm working a lot with semiconductor phyics, I wonder if there is a way to solve the common continuity equations together with the Poisson equation. Perhaps this is a known issue? I already tried it with pdepe, but continue to receive following error message:

Error using pdepe (line 293)
Spatial discretization has failed. Discretization supports only parabolic and elliptic equations, with flux term involving spatial
derivative.
Error in ContinuityEquations (line 27)
sol = pdepe(0,pdefun,icfun,bcfun,x,t);

I tried the following two versions of pdefun (all physical constants are set to 1). First version doesn't work because of the zero in the 3rd component of the flux-term I guess. Does anyone know what is wrong with the 2nd version?

function [c,f,s] = pdefun(x,t,u,DuDx)

    c = [1;1;0];

    f = [DuDx(1); DuDx(2); 0];

    s = [u(1)*DuDx(3)+u(3)*DuDx(1); ...
        -u(2)*DuDx(3)-u(3)*DuDx(2); ...
        -DuDx(3)+u(2)-u(1)];
end

function [c,f,s] = pdefun(x,t,u,DuDx)

    c = [1;1;0];

    f = [DuDx(1)-u(1)*DuDx(3); ...
        DuDx(2)+u(2)*DuDx(3); ...
        DuDx(3)];

    s = [-DuDx(3)*DuDx(1)+DuDx(1)*DuDx(3); ...
        DuDx(3)*DuDx(2)-DuDx(2)*DuDx(3); ...
        u(2)-u(1)];
end

Answer 1

[Sharmila Raghu]

The above error might occur if the boundary conditions are ill-posed. Please verify the boundary conditions to see if they are ill-posed. The boundary conditions specified as "p" and "q", follow this relationship:

p + q*f = 0

If "pr" and "qr" (the parameters for the right boundary) are both zero, this becomes 0+0=0, which is an ill-posed problem. To resolve this, please try setting "qr" to anything besides zero. This is equivalent to "f=0", which was probably the intended result.

Answer 2

[Sven B.]

Thanks for you answer! Indeed I had an ill-posed boundary condition, but unfortunately I still suffer from the same problem.

Find attached my complete code so far and the whole task in a better readable mathematical notation.

% Mesh

x = (-2:0.01:2)*1e-4; % cm
t = logspace(-16,-10,1000); % s

% Parameters

N_D = 1e17; % cm-3
N_A = 1e10; % cm-3
mu_n = 1000; % cm²/Vs
mu_p = 350; % cm²/Vs
epsilon_0 = 8.854187817e-14; % As/Vcm
epsilon_r = 11.9; 
e = 1.602177e-19; % As
k_B = 8.617330e-5; % eV/K

% Starting conditions

gauss = @(x,A,mu,sigma) A*exp(-(x-mu).^2/(2*sigma^2));

n0 = @(x) N_D + gauss(x,1e18,0,1e-5); % cm-3
p0 = @(x) zeros(1,numel(x))+N_A; % cm-3

dEdx0 = @(x) e/(epsilon_0*epsilon_r)*(p0(x)-n0(x)+N_D-N_A); % V/cm2
E0 = @(x) int(x,dEdx0)-0.5*integral(dEdx0,x(1),x(end)); % V/cm (see below for 'int')

% Solution

pdefun = @(x,t,u,DuDx) pdefunH(x,t,u,DuDx,N_D,N_A,mu_n,mu_p,epsilon_0,epsilon_r,e,k_B);
icfun = @(x) icfunH(x,n0,p0,E0);
bcfun = @(xl,ul,xr,ur,t) bcfunH(xl,ul,xr,ur,t,E0,x);

sol = pdepe(0,pdefun,icfun,bcfun,x,t);

function [c,f,s] = pdefunH(x,t,u,DuDx,N_D,N_A,mu_n,mu_p,epsilon_0,epsilon_r,e,k_B)

      c = [1;1;0];

      D_n = k_B*300/e*mu_n; %%cm2/s
      D_p = k_B*300/e*mu_p; %%cm2/s

      f = [D_n*DuDx(1); ...
          D_p*DuDx(2); ...
          eps*DuDx(3)];

      s = [mu_n*u(1)*DuDx(3)+mu_n*u(3)*DuDx(1); ...
          -mu_p*u(2)*DuDx(3)-mu_p*u(3)*DuDx(2); ...
          DuDx(3)-e/(epsilon_0*epsilon_r)*(u(2)-u(1)+N_D-N_A)];

end

function u = icfunH(x,n0,p0,E0)

      u = [n0(x);p0(x);E0(x)];

end

function [pl,ql,pr,qr] = bcfunH(xl,ul,xr,ur,t,E0,x)

      pl = [0;0;0];
      ql = [1;1;1];
      pr = [0;0;0];
      qr = [1;1;1];

end

function F = int( x, f )

      F = zeros(1,numel(x));

      for n = 1:numel(x)
         F(n) = integral(f,x(1),x(n)); 
      end

end