restart;
with(plots): Digits:= trunc(evalhf(Digits)); #generally a very efficient setting

# Here we solve a 1D problem in 3 regions. In each region, we have concentration and potential (c,phi) distributions,
# We first solve the unperturbed steady-state problem and then the linearized perturbation problem (which rely on the steady state).
# Each region is defined in x = 0..1, and the regions are connected by interface conditions that connect (c1(1),phi1(1)) to (c2(0),phi2(0))  and (c2(1),phi2(1)) to (c3(0),phi3(0))

Q:=10;   omega:=100;     J0:= 1.95;   # parameters

# The unperturbed steady-state

c1:=1-J0/2*x:               c3:=1-J0/2*(x-1):                   # concentration distributions in region 1 and 3    
c12:= eval(c1,x=1):        c32 := eval(c3,x=0):  
T1:=sqrt(Q^2+4*c12^2):     T3:=sqrt(Q^2+4*c32^2):           # the values of concentrations 1 and 3 at the interfaces with region 2
c21:=(T1-Q)/2:             c23:=(T3-Q)/2:                     # the values of concentration 2 at the interfaces with region 1 and 3 
I0:=fsolve(Q*i0/2/J0*ln((J0*T1-Q*i0)/(J0*T3-Q*i0))=(J0-T1+T3)/2,i0);   # the electrical current 
V:=(I0/J0+1)*ln(c32/c12)+ln((c21+Q)/(c23+Q))+(J0+2*c23-2*c21)/Q;     # the potential drop across the system 
c2:=solve(y-c21+Q*I0/2/J0*ln((Q*I0-Q*J0-2*J0*y)/(Q*I0-Q*J0-2*J0*c21))=-J0/2*x,y):  # concentration distribution in region 2 
phi1:=I0/J0*ln(c1)+V:   phi3:=I0/J0*ln(c3):                         # potential distribution in regions 1 and 3    
phi21:=I0/J0*ln(c12)+V-0.5*ln((c21+Q)/c21):    
phi2:=(2*c21-2*c2+Q*phi21-J0*x)/Q:      # potential distribution in region 2    

# The linearized problem 
# Unknowns: C11,C12,Phi11,Phi12,C21,C22,Phi21,Phi22,C31,C32,Phi31,Phi32,sigma1,sigma2 (sigma1 and sigma2 are constants along x)

#   Equations

# Region 1 Equations 

eq11 := omega*C11(x)-diff(diff(C12(x), x), x) = 0:                            
eqA1 := 2*c1*diff(Phi11(x), x)+2*(diff(phi1, x))*C11(x) = -sigma1: 
eq12 := omega*C12(x)+diff(diff(C11(x), x), x) = 0:                          
eqA2 := 2*c1*diff(Phi12(x), x)+2*(diff(phi1, x))*C12(x)=-sigma2:

 # Region 2 Equations 

eq21 := omega*C21(x)-diff(diff(C22(x), x)+Q/2*diff(Phi22(x), x), x)=0:      
eqB1 := 2*(c2+Q)*diff(Phi21(x), x)+2*(diff(phi2, x))*C21(x)=-sigma1:
eq22 :=  omega*C22(x)+diff(diff(C21(x), x)+Q/2*diff(Phi21(x), x), x) = 0:  
eqB2 := 2*(c2+Q)*diff(Phi22(x), x)+2*(diff(phi2, x))*C22(x)=-sigma2:

# Region 3 Equations 

eq31 := omega*C31(x)-diff(diff(C32(x), x), x)=0:                            
eqC1 := 2*c3*diff(Phi31(x), x)+2*(diff(phi3, x))*C31(x)=-sigma1:
eq32 := omega*C32(x)+diff(diff(C31(x), x), x) = 0:   
eqC2 := 2*c3*diff(Phi32(x), x)+2*(diff(phi3, x))*C32(x)=-sigma2:

EqSys := eq11, eq12, eq21, eq22, eq31, eq32, eqA1, eqA2, eqB1, eqB2, eqC1, eqC2;    # Equations system 

# Boundary conditions 

# Bcs at the outer ends of regions 1 and 3
Bc1 := C11(0) = 0, C12(0) = 0,  C31(1) = 0, C32(1) = 0, Phi11(0)=1, Phi12(0)=0, Phi31(1)=0, Phi32(1)=0:

# ECP continuity at the two interfaces (between region 1 and 2 and between 2 and 3) 
Intf1 := Phi21(0)-Phi11(1)=C11(1)/(eval(c1, x = 1))-C21(0)/(eval(c2, x = 0)+Q),
Phi22(0)-Phi12(1)=C12(1)/(eval(c1, x = 1))-C22(0)/(eval(c2, x = 0)+Q),
Phi21(0)-Phi11(1)=C21(0)/(eval(c2, x = 0))-C11(1)/(eval(c1, x = 1)),
Phi22(0)-Phi12(1)=C22(0)/(eval(c2, x = 0))-C12(1)/(eval(c1, x = 1)),
Phi21(1)-Phi31(0)=C31(0)/(eval(c3, x = 0))-C21(1)/(eval(c2, x = 1)+Q),
Phi22(1)-Phi32(0)=C32(0)/(eval(c3, x = 0))-C22(1)/(eval(c2, x = 1)+Q),
Phi21(1)-Phi31(0)=C21(1)/(eval(c2, x = 1))-C31(0)/(eval(c3, x = 0)),
Phi22(1)-Phi32(0)=C22(1)/(eval(c2, x = 1))-C32(0)/(eval(c3, x = 0)):

# Fluxes  continuity at the two interfaces (between region 1 and 2 and between 2 and 3)
Intf2 := (Q*sigma1+2*Q*D(phi2)(0)*C21(0))/(2*eval(c2, x = 0)+Q) = 2*D(C21)(0)-2*D(C11)(1),
(Q*sigma2+2*Q*D(phi2)(0)*C22(0))/(2*eval(c2, x = 0)+Q) = 2*D(C22)(0)-2*D(C12)(1),
(Q*sigma1+2*Q*D(phi2)(1)*C21(1))/(2*eval(c2, x = 1)+Q) = 2*D(C21)(1)-2*D(C31)(0),
(Q*sigma2+2*Q*D(phi2)(1)*C22(1))/(2*eval(c2, x = 1)+Q) = 2*D(C22)(1)-2*D(C32)(0): 

Bc := Bc1,Intf1,Intf2;

sys := {EqSys,Bc}:

sol1 := dsolve(sys, numeric, method=bvp[midrich],output = procedurelist);
Error, (in dsolve/numeric/bvp/convertsys) unable to convert to an explicit first-order system
Sigma1 := subs(sol1, sigma1);
Sigma2 := subs(sol1, sigma2);
Cond := Sigma1(0)+I*Sigma2(0);
ZR := Re(1/Cond);
ZI := Im(1/Cond);
X:=ZR,-ZI;