Why do I receive Error using pdepe (line 293) Spatial discretization has failed ?

I am trying to write a code to reproduce this model shown in the attached figure. C', A' and Phi' are u1, u2 and u3, respectively. Why do I receive Error using pdepe (line 293) Spatial discretization has failed ?

Here is my code.

clear, clc
global phio Z k T e B1 B2 B3 B4 B6 B7 B8 B9 B10 B11 
B1 = 5.55E-03 ; B2 = 1E-04 ; B3 = 8.33E-03 ; B4 = 1E-04 ;
B6 = 9.04E-011 ; B7 = B6; B8 = B4 ;  B9 = B3 ; B10 = B4 ;
B11 = B3 ; Z = 2 ; k = 1.38E-23 ; T = 300 ;  e = 1.6E-19 ; phio = 1 ;
function [c,f,s] = Nernest(x,t,u,dudx)
global phio Z k T e B1 B2 B3 B4 B6 B7 B8 B9 B10 B11 
c = [1; 1; 0];
f = [B1*u(1)*dudx(3)+B2*dudx(1) ; -B3*u(2)*dudx(3) + B4*dudx(2) ; -dudx(3) ] ;
s = [0; 0; B6*u(2)-B7*u(1)];
end
function u0 = Nernestic(x)
global phio Z k T e B1 B2 B3 B4 B6 B7 B8 B9 B10 B11 
u0 = [1;1;phio*x];
end
function [pl,ql,pr,qr] = Nernestbc(xl,ul,xr,ur,t)
global phio Z k T e B1 B2 B3 B4 B6 B7 B8 B9 B10 B11 
pl = [ 0      ;         0      ;  ul(3) + k*T*log(Z*ul(1))/(Z*e*phio)         ];
ql = [ 0      ;      0         ;   0                                        ];
pr = [ ur(1)  ;     ur(2)      ;  ur(3) - phio + (k*T*log(Z*ur(1))/(Z*e*phio))  ];        
qr = [  0     ;        0       ;  0                                         ];
end
x = linspace(0,1,25);
t = linspace(0,1,25);
m = 0;
sol = pdepe(m,@Nernest,@Nernestic,@Nernestbc,x,t);
u1 = sol(:,:,1);
u2 = sol(:,:,2);
u3 = sol(:,:,3);
plot(x,u3)

[Torsten]

The boundary conditions you set in your code are

0 = 0
0 = 0
phi(x=0)+ k*T*log(Z*C(x=0))/(Z*e*phio) = 0
C(x=1) = 0
A(x=1) = 0
phi(x=1)- phio + k*T*log(Z*C(x=1))/(Z*e*phio) = 0

Thus four of the six look wrong according to your mathematical description of the problem.

[Yaser]

Kindly, can you please help we with the correct code for the boundary conditions. I have no clue how to correct my mistakes in this code.

[Torsten]

Is C'(0) and C'(1) in your problem description the same as dC/dx @x=0 and dC/dx @x=1 ?

[Yaser]

In the problem the ' sign is just used to express dimensionless quantities not dervivatives. You can ignore any ' sign. I tried to hide them in the image to prevent any confusion.

The BCs are At x = 0

 du1/dx = 0

 -B8 * du2/dx + B9 * u2 * 
du3/dx =0

u3 = - kT ln (Z* u1)/ Z * e * phio

At x = 1

u1 = u2

-B10 * du2/dx + B11 * u2 * du3/dx =0

u3 = phio - kT ln (Z* u1)/ Z * e * phio

[Torsten]

The boundary conditions can include the value of the variable and/or the fluxes defined in the array f in your function "Nernest". The boundary conditions have the form

p(x,t,u) + f(x,t,u,du/dx)*q(x,t) = 0

The boundary condition for A, e.g., can directly be implemented in the code by setting

pl(2) = 0
ql(2) = 1
pr(2) = 0
qr(2) = 1

because this means

0 + 1*(-B3*ul(2)*dudx(3) + B4*dudx(2)) = 0 at x = 0
0 + 1*(-B3*ur(2)*dudx(3) + B4*dudx(2)) = 0 at x = 1.

thus exactly what you want to set.

The boundary condition for phi is set correctly in your code.

But I don't know how the boundary condition for C at x = 0, namely dC/dx = 0, could be set in the form

pl(1) + ql(1)*(B1*ul(1)*dudx(3)+B2*dudx(1)) = 0

with pl(1), ql(1) chosen suitably.

Thus it seems to me that the problem in its given form cannot be solved using "pdepe".

The elements marked by ?? are the problematic ones:

pl = [ ??            ;      0      ;  ul(3) + k*T*log(Z*ul(1))/(Z*e*phio)       ];
ql = [ ??            ;      1      ;  0                                        ];
pr = [ ur(1)-ur(2)   ;      0      ;  ur(3) - phio + (k*T*log(Z*ur(1))/(Z*e*phio))  ];        
qr = [  0            ;      1      ;  0                                         ];

[Yaser]

Thanks so much. I am so grateful to your answer 🙏.

[Torsten]

I made a slight mistake, but the answer remains the same (the problem seems unsolvable with "pdepe"):

But I don't know how the boundary condition for C at x = 0, namely dC/dx = 0, could be set in the form

pl(1) + ql(1)*(B1*ul(1)*dudx(3)+B2*dudx(1)) = dudx(1)

where pl(1), ql(1) had to be chosen suitably.