restart;
with(plots):
Digits:= trunc(evalhf(Digits));

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

# Steady-state solutions (unchanged)
c1:=1-J0/2*x:               
c3:=1-J0/2*(x-1):                   
c12:= eval(c1,x=1):        
c32 := eval(c3,x=0):  
T1:=sqrt(Q^2+4*c12^2):     
T3:=sqrt(Q^2+4*c32^2):           
c21:=(T1-Q)/2:             
c23:=(T3-Q)/2:                     
I0:=fsolve(Q*i0/2/J0*ln((J0*T1-Q*i0)/(J0*T3-Q*i0))=(J0-T1+T3)/2,i0);   
V:=(I0/J0+1)*ln(c32/c12)+ln((c21+Q)/(c23+Q))+(J0+2*c23-2*c21)/Q;     
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):                         
phi21:=I0/J0*ln(c12)+V-0.5*ln((c21+Q)/c21):    
phi2:=(2*c21-2*c2+Q*phi21-J0*x)/Q:      

# For the perturbation problem, let's define phi2 derivatives explicitly
phi2_deriv := diff(phi2, x);
phi2_deriv0 := eval(phi2_deriv, x=0);
phi2_deriv1 := eval(phi2_deriv, x=1);

# Interface concentration values
c1_1 := eval(c1, x=1);
c2_0 := eval(c2, x=0);
c2_1 := eval(c2, x=1);
c3_0 := eval(c3, x=0);

# Try a simplified approach - break into smaller systems
# Region 1 system
sys1 := {
    omega*C11(x) - diff(C12(x), x, x) = 0,
    omega*C12(x) + diff(C11(x), x, x) = 0,
    2*c1*diff(Phi11(x), x) + 2*diff(phi1, x)*C11(x) = -sigma1,
    2*c1*diff(Phi12(x), x) + 2*diff(phi1, x)*C12(x) = -sigma2,
    C11(0) = 0, C12(0) = 0, Phi11(0) = 1, Phi12(0) = 0
};

# Region 3 system  
sys3 := {
    omega*C31(x) - diff(C32(x), x, x) = 0,
    omega*C32(x) + diff(C31(x), x, x) = 0,
    2*c3*diff(Phi31(x), x) + 2*diff(phi3, x)*C31(x) = -sigma1,
    2*c3*diff(Phi32(x), x) + 2*diff(phi3, x)*C32(x) = -sigma2,
    C31(1) = 0, C32(1) = 0, Phi31(1) = 0, Phi32(1) = 0
};

# Alternative: Try solving the system without sigma1, sigma2 as unknowns first
# Convert to a first-order system manually

# For region 1, let's define new variables:
# Let U11 = diff(C11,x), U12 = diff(C12,x), V11 = diff(Phi11,x), V12 = diff(Phi12,x)
# Then we can write the system as 8 first-order ODEs for region 1

# However, given the complexity, you might need to:
# 1. Use a shooting method approach
# 2. Implement a finite difference method manually
# 3. Simplify the physics if possible

# For now, let's try a basic diagnostic:
print("System parameters defined successfully");
print("I0 =", I0);
print("V =", V);

# Check if c2 is properly defined
try
    c2_test := eval(c2, x=0.5);
    print("c2 at x=0.5:", c2_test);
catch:
    print("Error in c2 definition");
end try;