What does Error, (in dsolve/numeric/bvp) bad index into Matrix mean?
Also, I'm trying to run it, it is slow, any suggestions?

restart;
with(Student[VectorCalculus]);
with(DynamicSystems);
with(DEtools);
with(PDEtools, ReducedForm, declare, diff_table, dsubs);
NULL;
 #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:= 0.01;   # parameters
                            Q := 10

                          omega := 100

                           J0 := 0.01

# The unperturbed steady-state

c1:=1-J0/2*x: 
c3:=1-J0/2*(x-1):  
c12:= eval(c1,x=1); 
c32 := eval(c3,x=0); 
S1:=sqrt(Q^2+4*c12^2): 
S3:=sqrt(Q^2+4*c32^2):  
c21:=eval((S1-Q)/2); 
c23:=eval((S3-Q)/2);  
I0:=fsolve(Q*i0/2/J0*ln((J0*S1-Q*i0)/(J0*S3-Q*i0))=(J0-S1+S3)/2,i0);  
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):  
phi1:=I0/J0*ln(c1)+V: 
phi3:=I0/J0*ln(c3): 
dphi1:=diff(phi1,x); 
dphi3:=diff(phi3,x); 
phi21:=I0/J0*ln(c12)+V-0.5*ln((c21+Q)/c21); 
phi2:=(2*c21-2*c2+Q*phi21-J0*x)/Q: 
dphi2:=diff(phi2,x); 
dphi12 := eval(dphi1, x=1); 
dphi21 := eval(dphi2, x=0); 
dphi23 := eval(dphi2, x=1); 
dphi32 := eval(dphi3, x=0); 
INT1 := int(1/(2*c1), x = 0 .. 1); 
INT2 := int(1/(2*c2 + Q), x = 0 .. 1); 
INT3 := int(1/(2*c3), x = 0 .. 1); 
INT := INT1 + INT2 + INT3;
                      c12 := 0.9950000000

                       c32 := 1.005000000

                      c21 := 0.09804129000

                      c23 := 0.1000024500

                      I0 := 0.01419804328

                       V := 0.02539628566

                              0.007099021640   
                dphi1 := - --------------------
                           1 - 0.005000000000 x

                              0.007099021640        
           dphi3 := - ------------------------------
                      1.005000000 - 0.005000000000 x

                     phi21 := -2.299074561

dphi2 := (0.001000000000 LambertW(-0.2818670588 exp(-0.2818670588

   - 0.0007043224058 x)))/(1

   + LambertW(-0.2818670588 exp(-0.2818670588 - 0.0007043224058 x)

  )) - 0.001000000000

                   dphi12 := -0.007134695118

                   dphi21 := -0.001392499832

                   dphi23 := -0.001391964358

                   dphi32 := -0.007063703124

                      INT1 := 0.5012541824

                     INT2 := 0.09805801917

                      INT3 := 0.4987541511

                       INT := 1.098066353

sys1 := {
-omega*C11(x)+diff(diff(C12(x), x), x)=0,
omega*C12(x)+diff(diff(C11(x), x), x) = 0,
-omega*C21(x)+diff(diff(C22(x), x)+(c2*sigma2-C22(x)*dphi2*Q)/(2*c2+Q), x) =0,
 omega*C22(x)+diff(diff(C21(x), x)+(c2*sigma1-C21(x)*dphi2*Q)/(2*c2+Q), x) = 0,
-omega*C31(x)+diff(diff(C32(x), x), x)=0,
omega*C32(x)+diff(diff(C31(x), x), x) = 0
}:

sys2 := {
diff(FA1(x), x) = C11(x)*dphi1/c1,
diff(FA2(x), x) = C21(x)*dphi2/(c2+Q/2),
diff(FA3(x), x) = C31(x)*dphi3/c3,
diff(FB1(x), x) = C12(x)*dphi1/c1,
diff(FB2(x), x) = C22(x)*dphi2/(c2+Q/2),
diff(FB3(x), x) = C32(x)*dphi3/c3
}: 

Bc := {
C11(0) = 0, C12(0) = 0,  C31(1) = 0, C32(1) = 0,
FA1(0) = 0, FB1(0) = 0,  FA3(1) = 0, FB3(1) = 0, 

2*C11(1)/c12 = C21(0)/(c21+Q)+C21(0)/c21, 
2*C12(1)/c12 = C22(0)/(c21+Q)+C22(0)/c21,
C21(1)/(c23+Q)+C21(1)/c23 = 2*C31(0)/c32,
C22(1)/(c23+Q)+C22(1)/c23 = 2*C32(0)/c32,

D(C11)(1)+dphi12*C11(1)-sigma1/2-c12*D(FA1)(1) = D(C21)(0)+dphi21*C21(0)-(c21+Q)*sigma1/(2*c21+Q)-(c21+Q)*D(FA2)(0),
D(C12)(1)+dphi12*C12(1)-sigma2/2-c12*D(FB1)(1) = D(C22)(0)+dphi21*C22(0)-(c21+Q)*sigma2/(2*c21+Q)-(c21+Q)*D(FB2)(0),
D(C11)(1)-dphi12*C11(1)+sigma1/2+c12*D(FA1)(1) = D(C21)(0)-dphi21*C21(0)+c21*sigma1/(2*c21+Q)+c21*D(FA2)(0),
D(C12)(1)-dphi12*C12(1)+sigma2/2+c12*D(FB1)(1) = D(C22)(0)-dphi21*C22(0)+c21*sigma2/(2*c21+Q)+c21*D(FB2)(0),

D(C31)(0)+dphi32*C31(0)-sigma1/2-c32*D(FA3)(0) = D(C21)(1)+dphi23*C21(1)-(c23+Q)*sigma1/(2*c23+Q)-(c23+Q)*D(FA2)(1),
D(C32)(0)+dphi32*C32(0)-sigma2/2-c32*D(FB3)(0) = D(C22)(1)+dphi23*C22(1)-(c23+Q)*sigma2/(2*c23+Q)-(c23+Q)*D(FB2)(1),
D(C31)(0)-dphi32*C31(0)+sigma1/2+c32*D(FA3)(0) = D(C21)(1)-dphi23*C21(1)+c23*sigma1/(2*c23+Q)+c23*D(FA2)(1),
D(C32)(0)-dphi32*C32(0)+sigma2/2+c32*D(FB3)(0) = D(C22)(1)-dphi23*C22(1)+c23*sigma2/(2*c23+Q)+c23*D(FB2)(1)
}:
 
 

all_sys := sys1 union sys2 union Bc:
sol1 := dsolve(all_sys, initmesh = 100, maxmesh = 15000, numeric, method = bvp[midrich], output = listprocedure):
#(all_sys, numeric, method = bvp[midrich]);

Error, (in dsolve/numeric/bvp) bad index into Matrix

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;