# Equation that governs the steady front positions
F1:=1/(R*S0+1-S0):
F:=F1*(1-R*S0*R3*k*(1-R1)/(C*Pec_i)):
F2:=Pec_i*F:
frontlocationeq:= Theta0/C=(exp(Pec_i*F*(1-S0))-1)*(H+1/(1-exp(-C/k*Pec_i*F*S0))):

# Dispersion equation
N0:=F1*C/k*Pec_i/R*((R-1)+k*R*R3*(1-R1)/(C*Pec_i)):
omega0:=1/k*Pec_i*F*C:
Omega0:=Pec_i*F:
N1:=omega0/(exp(omega0*S0)-1):
N3:=N1*exp(omega0*S0):
N4:=N2*exp(Omega0*(S0-1)):
N2:=Omega0/(exp(Omega0*(S0-1))-1):
C1 := 1/2*((-k*R*R3*(-1+R1)+C*Pec_i*(-1+R))*F1*l*cosh(l*(-1+S0))-k*sinh(l*(-1+S0))*variphi*sigma*(-1+R1))*Phi/k/l/(cosh(l*(-1+S0))*sinh(l*S0)*R-sinh(l*(-1+S0))*cosh(l*S0)): C2:=1/2*Phi*(F1*l*((-R3*(-1+R1)*k+C*Pec_i)*R-C*Pec_i)*cosh(l*S0)-R*k*sinh(l*S0)*sigma*variphi*(-1+R1))*(exp(-l))/R/l/k/(cosh(l*(-1+S0))*sinh(l*S0)*R-sinh(l*(-1+S0))*cosh(l*S0)):
phi2(x):=-4*(exp(1/2*Omega0*x)*(k*C2*R*exp(1/2*Omega0*S0)*l*N2*(exp(-Omega0)*(C*Omega0*l+k*R1*E2*sigma)+(-C*Omega0*l+k*R1*E2*sigma)*exp(-Omega0+2*l))*cosh(1/2*S0*((Omega0^2+4*l^2)*C^2+4*C*k*R1*E2*sigma)^(1/2)/C)+k*R*l*C2*N2*(C*Omega0*l-k*R1*E2*sigma)*exp((S0-1)*Omega0-l*(-2+S0))+(C*Omega0*l+k*R1*E2*sigma)*(-C2*N2*R*k*exp((S0-1)*Omega0+l*S0)*l+Phi*N4*(C*Omega0*l-k*R1*E2*sigma)))*sinh(1/2*x*((Omega0^2+4*l^2)*C^2+4*C*k*R1*E2*sigma)^(1/2)/C)-k*C2*R*sinh(1/2*S0*((Omega0^2+4*l^2)*C^2+4*C*k*R1*E2*sigma)^(1/2)/C)*exp(1/2*Omega0*S0)*(exp(1/2*Omega0*x)*(exp(-Omega0)*(C*Omega0*l+k*R1*E2*sigma)+(-C*Omega0*l+k*R1*E2*sigma)*exp(-Omega0+2*l))*cosh(1/2*x*((Omega0^2+4*l^2)*C^2+4*C*k*R1*E2*sigma)^(1/2)/C)+(C*Omega0*l-k*R1*E2*sigma)*exp(Omega0*(x-1)-l*(-2+x))-(C*Omega0*l+k*R1*E2*sigma)*exp(Omega0*(x-1)+l*x))*l*N2)*(k*R1*E2*sigma+C*(1/4*Omega0^2+l^2))/(C*Omega0^2+4*k*R1*E2*sigma+4*l^2*C)/sinh(1/2*S0*((Omega0^2+4*l^2)*C^2+4*C*k*R1*E2*sigma)^(1/2)/C)/exp(1/2*Omega0*S0)/(C^2*Omega0^2*l^2-E2^2*R1^2*k^2*sigma^2):
phi1(x) := (omega0*l-E1*sigma)*N1*(exp(1/2*omega0*x)*(-2*l*C1*E1*sigma*exp(1/2*omega0*S0)*cosh(1/2*S0*(omega0^2+4*E1*sigma+4*l^2)^(1/2))-l*C1*(omega0*l-E1*sigma)*exp(-(-omega0+l)*S0)+(l*C1*exp((omega0+l)*S0)-Phi*exp(omega0*S0)*(omega0*l-E1*sigma))*(omega0*l+E1*sigma))*sinh(1/2*x*(omega0^2+4*E1*sigma+4*l^2)^(1/2))+(2*E1*sigma*exp(1/2*omega0*x)*cosh(1/2*x*(omega0^2+4*E1*sigma+4*l^2)^(1/2))+(omega0*l-E1*sigma)*exp((omega0-l)*x)-exp((omega0+l)*x)*(omega0*l+E1*sigma))*sinh(1/2*S0*(omega0^2+4*E1*sigma+4*l^2)^(1/2))*l*C1*exp(1/2*omega0*S0))*(omega0*l+E1*sigma)/(l^2*omega0^2-E1^2*sigma^2)^2/sinh(1/2*S0*(omega0^2+4*E1*sigma+4*l^2)^(1/2))/exp(1/2*omega0*S0): 
# In A, cosh(l)-sinh(l) replaced with exp(-l) wherever it appears
A:=variphi*sigma*H*Phi-((-Phi*(omega0*l-E1*sigma)*(-1/4*exp(1/2*Omega0*S0)*(F1*(-k*R*R3*(-1+R1)+C*Pec_i*(-1+R))*l*cosh(l*S0)-R*k*sinh(l*S0)*sigma*variphi*(-1+R1))*Theta0*N2*(exp(-l))*(omega0*l-E1*sigma)*((Omega0^2+4*l^2)*C^2+4*C*k*R1*E2*sigma)^(1/2)*sinh(1/2*S0*(omega0^2+4*E1*sigma+4*l^2)^(1/2))*(exp(-Omega0)*(C*Omega0*l+k*R1*E2*sigma)+(-C*Omega0*l+k*R1*E2*sigma)*exp(-Omega0+2*l))*(omega0*l+E1*sigma)*sinh(1/2*S0*((Omega0^2+4*l^2)*C^2+4*C*k*R1*E2*sigma)^(1/2)/C)^2+C*(-1/2*sigma*N1*((-k*R*R3*(-1+R1)+C*Pec_i*(-1+R))*F1*l*cosh(l*(-1+S0))-k*sinh(l*(-1+S0))*variphi*sigma*(-1+R1))*(omega0^2+4*E1*sigma+4*l^2)^(1/2)*E1*(C*Omega0*l-k*R1*E2*sigma)*exp(1/2*omega0*S0)*(C*Omega0*l+k*R1*E2*sigma)*sinh(1/2*S0*(omega0^2+4*E1*sigma+4*l^2)^(1/2))^2+(1/2*(F1*(-k*R*R3*(-1+R1)+C*Pec_i*(-1+R))*l*cosh(l*S0)-R*k*sinh(l*S0)*sigma*variphi*(-1+R1))*Theta0*N2*(exp(-l))*(omega0*l-E1*sigma)*(-1/2*Omega0+l)*(C*Omega0*l-k*R1*E2*sigma)*(omega0*l+E1*sigma)*exp((-1+S0)*Omega0-l*(-2+S0))+(C*Omega0*l+k*R1*E2*sigma)*(1/2*(F1*(-k*R*R3*(-1+R1)+C*Pec_i*(-1+R))*l*cosh(l*S0)-R*k*sinh(l*S0)*sigma*variphi*(-1+R1))*Theta0*N2*(exp(-l))*(omega0*l-E1*sigma)*(l+1/2*Omega0)*(omega0*l+E1*sigma)*exp((Omega0+l)*S0-Omega0)+(1/4*N1*omega0*((-k*R*R3*(-1+R1)+C*Pec_i*(-1+R))*F1*l*cosh(l*(-1+S0))-k*sinh(l*(-1+S0))*variphi*sigma*(-1+R1))*(omega0*l-E1*sigma)*exp(-(-omega0+l)*S0)+1/2*N1*((-k*R*R3*(-1+R1)+C*Pec_i*(-1+R))*F1*l*cosh(l*(-1+S0))-k*sinh(l*(-1+S0))*variphi*sigma*(-1+R1))*(omega0*l-E1*sigma)*(-omega0+l)*exp((omega0-l)*S0)+(omega0*l+E1*sigma)*((1/2*F1*(-k*R*R3*(-1+R1)+C*Pec_i*(-1+R))*N1*(l+1/2*omega0)*l*exp((omega0+l)*S0)+(F1*(-k*R*R3*(-1+R1)+C*Pec_i*(-1+R))*H*l*cosh(l*S0)+1/2*R*sinh(l*S0)*(omega0*N1*exp(omega0*S0)*k-2*omega0*N3*k-Omega0*Theta0*N4))*(omega0*l-E1*sigma))*cosh(l*(-1+S0))-(1/2*sigma*N1*variphi*(-1+R1)*(l+1/2*omega0)*k*exp((omega0+l)*S0)+(omega0*l-E1*sigma)*cosh(l*S0)*(1/2*omega0*N1*exp(omega0*S0)*k+H*k*variphi*(-1+R1)*sigma-omega0*N3*k-1/2*Omega0*Theta0*N4))*sinh(l*(-1+S0))))*(C*Omega0*l-k*R1*E2*sigma)))*sinh(1/2*S0*(omega0^2+4*E1*sigma+4*l^2)^(1/2))-1/4*(-2*sigma*exp(1/2*omega0*S0)*E1*((-k*R*R3*(-1+R1)+C*Pec_i*(-1+R))*F1*l*cosh(l*(-1+S0))-k*sinh(l*(-1+S0))*variphi*sigma*(-1+R1))*cosh(1/2*S0*(omega0^2+4*E1*sigma+4*l^2)^(1/2))-(omega0*l-E1*sigma)*((-k*R*R3*(-1+R1)+C*Pec_i*(-1+R))*F1*l*cosh(l*(-1+S0))-k*sinh(l*(-1+S0))*variphi*sigma*(-1+R1))*exp(-(-omega0+l)*S0)+(((-k*R*R3*(-1+R1)+C*Pec_i*(-1+R))*l*F1*exp((omega0+l)*S0)-2*R*k*exp(omega0*S0)*sinh(l*S0)*(omega0*l-E1*sigma))*cosh(l*(-1+S0))+2*sinh(l*(-1+S0))*(-1/2*variphi*sigma*(-1+R1)*exp((omega0+l)*S0)+exp(omega0*S0)*cosh(l*S0)*(omega0*l-E1*sigma))*k)*(omega0*l+E1*sigma))*N1*(omega0^2+4*E1*sigma+4*l^2)^(1/2)*(C*Omega0*l-k*R1*E2*sigma)*cosh(1/2*S0*(omega0^2+4*E1*sigma+4*l^2)^(1/2))*(C*Omega0*l+k*R1*E2*sigma))*sinh(1/2*S0*((Omega0^2+4*l^2)*C^2+4*C*k*R1*E2*sigma)^(1/2)/C)+1/4*Theta0*(omega0*l-E1*sigma)*((Omega0^2+4*l^2)*C^2+4*C*k*R1*E2*sigma)^(1/2)*sinh(1/2*S0*(omega0^2+4*E1*sigma+4*l^2)^(1/2))*(exp(1/2*Omega0*S0)*(F1*(-k*R*R3*(-1+R1)+C*Pec_i*(-1+R))*l*cosh(l*S0)-R*k*sinh(l*S0)*sigma*variphi*(-1+R1))*N2*(exp(-l))*(exp(-Omega0)*(C*Omega0*l+k*R1*E2*sigma)+(-C*Omega0*l+k*R1*E2*sigma)*exp(-Omega0+2*l))*cosh(1/2*S0*((Omega0^2+4*l^2)*C^2+4*C*k*R1*E2*sigma)^(1/2)/C)+(F1*(-k*R*R3*(-1+R1)+C*Pec_i*(-1+R))*l*cosh(l*S0)-R*k*sinh(l*S0)*sigma*variphi*(-1+R1))*N2*(exp(-l))*(C*Omega0*l-k*R1*E2*sigma)*exp((-1+S0)*Omega0-l*(-2+S0))-((F1*(-k*R*R3*(-1+R1)+C*Pec_i*(-1+R))*l*cosh(l*S0)-R*k*sinh(l*S0)*sigma*variphi*(-1+R1))*N2*(exp(-l))*exp((Omega0+l)*S0-Omega0)-2*N4*(cosh(l*(-1+S0))*sinh(l*S0)*R-sinh(l*(-1+S0))*cosh(l*S0))*(C*Omega0*l-k*R1*E2*sigma))*(C*Omega0*l+k*R1*E2*sigma))*(omega0*l+E1*sigma)*cosh(1/2*S0*((Omega0^2+4*l^2)*C^2+4*C*k*R1*E2*sigma)^(1/2)/C))*(omega0*l+E1*sigma)/k/(cosh(l*(-1+S0))*sinh(l*S0)*R-sinh(l*(-1+S0))*cosh(l*S0))/(l^2*omega0^2-E1^2*sigma^2)^2/sinh(1/2*S0*(omega0^2+4*E1*sigma+4*l^2)^(1/2))/C/(-R1^2*E2^2*sigma^2*k^2+C^2*Omega0^2*l^2)/sinh(1/2*S0*((Omega0^2+4*l^2)*C^2+4*C*k*R1*E2*sigma)^(1/2)/C))):