## 220 Reputation

17 years, 28 days
Glasgow, United Kingdom

## Re;...

I want to get the data fromt the plot in a .txt file in the form of two columns.

## Re...

I agree with you. Its quite clear that in the case of 2 and 3 F(S) is always positive. Actually,  its the question of 1 and 4. What I want is to present the whole nature of F(S) in a plot form for all these cases. That plot should be very clearly mentioning the nature of F(S) for those 4 cases.

Thanks

## Re...

I agree with you. Its quite clear that in the case of 2 and 3 F(S) is always positive. Actually,  its the question of 1 and 4. What I want is to present the whole nature of F(S) in a plot form for all these cases. That plot should be very clearly mentioning the nature of F(S) for those 4 cases.

Thanks

## re...

still waiting for any response.

thx

## RE: Complex plot....

How to plot a set of complex data points.

## RE: RootFinding[Analytic]...

I tried but its not working. Well, Lets wait may be some one come up with a better idea. Anyway thanks

## RE: RootFinding[Analytic]...

I tried but its not working. Well, Lets wait may be some one come up with a better idea. Anyway thanks

## Re...

 > restart:with(MultiSeries):
 > Digits:=10:
 > variphi:=0.38:R1:=0.0009:R2:=8.75:E1:=1:E2:=1:R4:=177.6:C:=1.96:k:=4:H:=7:Pec_i:=-120:Phi:=1:Q:=1.9:R:=20:
 > A:=(Q*(R*S+1-S)*(H+1/(1-exp(Pec_i)))*exp(C*Pec_i*S/(k*(R*S+1-S)))-H)*(1-exp(Pec_i*(1-S)/(R*S+1-S)))=1;
 (1)
 > eqn2:=Student:-Calculus1:-Roots(subs(H=7,R=20,C=1.96,k=4,Pec_i=-120,Q=1.9,A),S);
 (2)
 > S0:=.3592598898850950896270594340253413222287;
 (3)
 > F:=1/(R*S0+1-S0);
 (4)
 > w_0:=1/k*Pec_i*F*C;
 (5)
 > w_1:=Pec_i*F;
 (6)
 > M_1:=evalf(Q*Pec_i/k*(H+1/(1-exp(Pec_i))));
 (7)
 > M_2:=M_1*exp(w_0*S0);
 (8)
 > M_3:=w_1/(1-exp(w_1*(S0-1)));
 (9)
 > M_4:=M_3*exp(w_1*(S0-1));
 (10)

R3:=-(-cosh(1/2*S0*(w_0^2)^(1/2))*(w_0^2)^(1/2)+sinh(1/2*S0*(w_0^2)^(1/2))*w_0)*(S0*R2*R1+1-S0)*w_1*((M_2*k*w_0+1/2*M_4*w_1)*C*sinh(1/2*(-1+S0)*(C^2*w_1^2)^(1/2)/C)-1/2*(C^2*w_1^2)^(1/2)*cosh(1/2*(-1+S0)*(C^2*w_1^2)^(1/2)/C)*M_4)*w_0/((-(1/2*M_3*C*(2+w_1*(-1+S0))*w_0*exp(w_1*(-1+S0))+R1*(C*M_1*(S0*w_0+1)*exp(S0*w_0)-w_0*H)*k*w_1)*(w_0^2)^(1/2)*cosh(1/2*S0*(w_0^2)^(1/2))+(1/2*M_3*C*(2+w_1*(-1+S0))*w_0*exp(w_1*(-1+S0))+R1*(C*M_1*(S0*w_0+1)*exp(S0*w_0)-w_0*H)*k*w_1)*w_0*sinh(1/2*S0*(w_0^2)^(1/2))+(w_0^2)^(1/2)*k*M_1*exp(1/2*S0*w_0)*w_1*C*R1)*sinh(1/2*(-1+S0)*(C^2*w_1^2)^(1/2)/C)-1/2*M_3*w_0*(C^2*w_1^2)^(1/2)*cosh(1/2*(-1+S0)*(C^2*w_1^2)^(1/2)/C)*exp(w_1*(-1+S0))*(-cosh(1/2*S0*(w_0^2)^(1/2))*(w_0^2)^(1/2)+sinh(1/2*S0*(w_0^2)^(1/2))*w_0)*(-1+S0))/(R1-1)/R2;

 (11)
 > R3:=-38.16;
 (12)
 > A2:=variphi*sigma*Phi = -2*(((1/4*(cosh(l)-sinh(l))*(R3*cosh(l*S0)*l+sinh(l*S0)*variphi*sigma)*w_0*(-l*w_0+E1*sigma)*(R1-1)*C^2*(l*C*w_1-sigma*E2*R1*k)*(l*w_0+E1*sigma)*M_3*R2*(-1/2*w_1+l)*exp(S0*(w_1-l)-w_1+2*l)+(1/4*(1/2*w_1+l)*(cosh(l)-sinh(l))*(R3*cosh(l*S0)*l+sinh(l*S0)*variphi*sigma)*w_0*(-l*w_0+E1*sigma)*(R1-1)*C^2*(l*w_0+E1*sigma)*M_3*R2*exp(S0*(w_1+l)-w_1)+(-1/2*(l^2-1/2*l*w_0+E1*sigma+1/2*w_0^2)*(-l*w_0+E1*sigma)*(R1-1)*(l*R3*R2*cosh(l*(-1+S0))*R1+variphi*sinh(l*(-1+S0))*sigma)*C*M_1*k*exp(-(-w_0+l)*S0)+((-1/2*l*(R1-1)*C*R3*M_1*(l^2+1/2*l*w_0+E1*sigma+1/2*w_0^2)*k*exp((w_0+l)*S0)+(-l*w_0+E1*sigma)*(-1/2*H*k*R3*l*w_0*(R1-1)*cosh(l*S0)+(M_2*k*sigma*E1+1/2*M_2*w_0^2*k+l^2*M_2*k+1/4*w_1*M_4*w_0)*C*sinh(l*S0)))*R1*R2*cosh(l*(-1+S0))-(1/2*(R1-1)*sigma*variphi*C*M_1*(l^2+1/2*l*w_0+E1*sigma+1/2*w_0^2)*k*exp((w_0+l)*S0)+((M_2*k*sigma*E1+1/2*M_2*w_0^2*k+l^2*M_2*k+1/4*w_1*M_4*w_0)*C+1/2*H*k*sigma*variphi*w_0*(R1-1))*(-l*w_0+E1*sigma)*cosh(l*S0))*sinh(l*(-1+S0)))*(l*w_0+E1*sigma))*(l*C*w_1-sigma*E2*R1*k))*(l*C*w_1+sigma*E2*R1*k))*sinh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))-1/2*(w_0^2+4*E1*sigma+4*l^2)^(1/2)*((1/2*(cosh(l)-sinh(l))*(R3*cosh(l*S0)*l+sinh(l*S0)*variphi*sigma)*(-l*w_0+E1*sigma)*(R1-1)*C^2*(l*C*w_1-sigma*E2*R1*k)*(l*w_0+E1*sigma)*M_3*R2*(-1/2*w_1+l)*exp(S0*(w_1-l)-w_1+2*l)+(l*C*w_1+sigma*E2*R1*k)*(1/2*(1/2*w_1+l)*(cosh(l)-sinh(l))*(R3*cosh(l*S0)*l+sinh(l*S0)*variphi*sigma)*(-l*w_0+E1*sigma)*(R1-1)*C^2*(l*w_0+E1*sigma)*M_3*R2*exp(S0*(w_1+l)-w_1)+(1/2*C*k*M_1*(-w_0+l)*(R1-1)*(l*R3*R2*cosh(l*(-1+S0))*R1+variphi*sinh(l*(-1+S0))*sigma)*(-l*w_0+E1*sigma)*exp(-(-w_0+l)*S0)+((-1/2*C*k*R3*M_1*l*(w_0+l)*(R1-1)*exp((w_0+l)*S0)+(-l*w_0+E1*sigma)*(-H*k*R3*l*(R1-1)*cosh(l*S0)+C*(M_2*k*w_0+1/2*M_4*w_1)*sinh(l*S0)))*R1*R2*cosh(l*(-1+S0))-sinh(l*(-1+S0))*(1/2*C*k*M_1*sigma*variphi*(w_0+l)*(R1-1)*exp((w_0+l)*S0)+(-l*w_0+E1*sigma)*cosh(l*S0)*((M_2*k*w_0+1/2*M_4*w_1)*C+H*k*sigma*variphi*(R1-1))))*(l*w_0+E1*sigma))*(l*C*w_1-sigma*E2*R1*k)))*cosh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))+C*k*M_1*w_0*exp(1/2*S0*w_0)*(R1-1)*(l*C*w_1-sigma*E2*R1*k)*(l*C*w_1+sigma*E2*R1*k)*(l*R3*R2*cosh(l*(-1+S0))*R1+variphi*sinh(l*(-1+S0))*sigma)*(E1*sigma+l^2)))*sinh(1/2*(-1+S0)*((w_1^2+4*l^2)*C^2+4*C*sigma*E2*R1*k)^(1/2)/C)-1/8*(-cosh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))*(w_0^2+4*E1*sigma+4*l^2)^(1/2)+sinh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))*w_0)*((w_1^2+4*l^2)*C^2+4*C*sigma*E2*R1*k)^(1/2)*(-l*w_0+E1*sigma)*(l*w_0+E1*sigma)*((-C*M_3*R2*(cosh(l)-sinh(l))*(R1-1)*(R3*cosh(l*S0)*l+sinh(l*S0)*variphi*sigma)*(l*C*w_1-sigma*E2*R1*k)*exp(S0*(w_1-l)-w_1+2*l)+(C*M_3*R2*(cosh(l)-sinh(l))*(R1-1)*(R3*cosh(l*S0)*l+sinh(l*S0)*variphi*sigma)*exp(S0*(w_1+l)-w_1)+2*M_4*(R1*R2*cosh(l*(-1+S0))*sinh(l*S0)-cosh(l*S0)*sinh(l*(-1+S0)))*(l*C*w_1-sigma*E2*R1*k))*(l*C*w_1+sigma*E2*R1*k))*cosh(1/2*(-1+S0)*((w_1^2+4*l^2)*C^2+4*C*sigma*E2*R1*k)^(1/2)/C)-2*C*k*E2*M_3*R1*sigma*R2*exp(1/2*w_1*(-1+S0)+l)*(cosh(l)-sinh(l))*(R1-1)*(R3*cosh(l*S0)*l+sinh(l*S0)*variphi*sigma)))*Phi/H/k/(l^2*C^2*w_1^2-sigma^2*E2^2*R1^2*k^2)/sinh(1/2*(-1+S0)*((w_1^2+4*l^2)*C^2+4*C*sigma*E2*R1*k)^(1/2)/C)/(R1*R2*cosh(l*(-1+S0))*sinh(l*S0)-cosh(l*S0)*sinh(l*(-1+S0)))/(-cosh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))*(w_0^2+4*E1*sigma+4*l^2)^(1/2)+sinh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))*w_0)/(-l^2*w_0^2+E1^2*sigma^2):
 > fsolve({diff((rhs-lhs)(A2),sigma)=0,(rhs-lhs)(A2)=0},{sigma,l=0.5});
 (13)
 > plots:-animate(plot,[[rhs,lhs](A2),sigma=-350..5,-130..10],l=1.5..3.5,frames=20);
 > plots:-animate(plot,[[rhs,lhs](A2),sigma=-350..5,-130..10],l=0..4.5,frames=20);
 > fsolve( [subs(l=4,A2)],sigma=-50, complex);
 >
 (14)
 > fsolve( [subs(l=3.9,A2)],sigma=-51.34547980+1e-5*I, complex);
 (15)
 >
 > fsolve( [subs(l=3.64,A2)],sigma=-22.44118854+1e-5*I, complex);
 (16)
 > eqn2:=Student:-Calculus1:-Roots(subs(l=3.5,A2),sigma);
 (17)
 >
 >
 >
 >

Thanks man. Thats what I was looking for.

But my friend, there is still a problem. In the interval [3.5, 8], sigma has multi real soultions. But the in interval [0, 3.5],  there is a big jump. Now i am looking to find that,  Is there any complex solutions for sigma for `l=0 to 3.5'.

## Re...

 > restart:with(MultiSeries):
 > Digits:=10:
 > variphi:=0.38:R1:=0.0009:R2:=8.75:E1:=1:E2:=1:R4:=177.6:C:=1.96:k:=4:H:=7:Pec_i:=-120:Phi:=1:Q:=1.9:R:=20:
 > A:=(Q*(R*S+1-S)*(H+1/(1-exp(Pec_i)))*exp(C*Pec_i*S/(k*(R*S+1-S)))-H)*(1-exp(Pec_i*(1-S)/(R*S+1-S)))=1;
 (1)
 > eqn2:=Student:-Calculus1:-Roots(subs(H=7,R=20,C=1.96,k=4,Pec_i=-120,Q=1.9,A),S);
 (2)
 > S0:=.3592598898850950896270594340253413222287;
 (3)
 > F:=1/(R*S0+1-S0);
 (4)
 > w_0:=1/k*Pec_i*F*C;
 (5)
 > w_1:=Pec_i*F;
 (6)
 > M_1:=evalf(Q*Pec_i/k*(H+1/(1-exp(Pec_i))));
 (7)
 > M_2:=M_1*exp(w_0*S0);
 (8)
 > M_3:=w_1/(1-exp(w_1*(S0-1)));
 (9)
 > M_4:=M_3*exp(w_1*(S0-1));
 (10)

R3:=-(-cosh(1/2*S0*(w_0^2)^(1/2))*(w_0^2)^(1/2)+sinh(1/2*S0*(w_0^2)^(1/2))*w_0)*(S0*R2*R1+1-S0)*w_1*((M_2*k*w_0+1/2*M_4*w_1)*C*sinh(1/2*(-1+S0)*(C^2*w_1^2)^(1/2)/C)-1/2*(C^2*w_1^2)^(1/2)*cosh(1/2*(-1+S0)*(C^2*w_1^2)^(1/2)/C)*M_4)*w_0/((-(1/2*M_3*C*(2+w_1*(-1+S0))*w_0*exp(w_1*(-1+S0))+R1*(C*M_1*(S0*w_0+1)*exp(S0*w_0)-w_0*H)*k*w_1)*(w_0^2)^(1/2)*cosh(1/2*S0*(w_0^2)^(1/2))+(1/2*M_3*C*(2+w_1*(-1+S0))*w_0*exp(w_1*(-1+S0))+R1*(C*M_1*(S0*w_0+1)*exp(S0*w_0)-w_0*H)*k*w_1)*w_0*sinh(1/2*S0*(w_0^2)^(1/2))+(w_0^2)^(1/2)*k*M_1*exp(1/2*S0*w_0)*w_1*C*R1)*sinh(1/2*(-1+S0)*(C^2*w_1^2)^(1/2)/C)-1/2*M_3*w_0*(C^2*w_1^2)^(1/2)*cosh(1/2*(-1+S0)*(C^2*w_1^2)^(1/2)/C)*exp(w_1*(-1+S0))*(-cosh(1/2*S0*(w_0^2)^(1/2))*(w_0^2)^(1/2)+sinh(1/2*S0*(w_0^2)^(1/2))*w_0)*(-1+S0))/(R1-1)/R2;

 (11)
 > R3:=-38.16;
 (12)
 > A2:=variphi*sigma*Phi = -2*(((1/4*(cosh(l)-sinh(l))*(R3*cosh(l*S0)*l+sinh(l*S0)*variphi*sigma)*w_0*(-l*w_0+E1*sigma)*(R1-1)*C^2*(l*C*w_1-sigma*E2*R1*k)*(l*w_0+E1*sigma)*M_3*R2*(-1/2*w_1+l)*exp(S0*(w_1-l)-w_1+2*l)+(1/4*(1/2*w_1+l)*(cosh(l)-sinh(l))*(R3*cosh(l*S0)*l+sinh(l*S0)*variphi*sigma)*w_0*(-l*w_0+E1*sigma)*(R1-1)*C^2*(l*w_0+E1*sigma)*M_3*R2*exp(S0*(w_1+l)-w_1)+(-1/2*(l^2-1/2*l*w_0+E1*sigma+1/2*w_0^2)*(-l*w_0+E1*sigma)*(R1-1)*(l*R3*R2*cosh(l*(-1+S0))*R1+variphi*sinh(l*(-1+S0))*sigma)*C*M_1*k*exp(-(-w_0+l)*S0)+((-1/2*l*(R1-1)*C*R3*M_1*(l^2+1/2*l*w_0+E1*sigma+1/2*w_0^2)*k*exp((w_0+l)*S0)+(-l*w_0+E1*sigma)*(-1/2*H*k*R3*l*w_0*(R1-1)*cosh(l*S0)+(M_2*k*sigma*E1+1/2*M_2*w_0^2*k+l^2*M_2*k+1/4*w_1*M_4*w_0)*C*sinh(l*S0)))*R1*R2*cosh(l*(-1+S0))-(1/2*(R1-1)*sigma*variphi*C*M_1*(l^2+1/2*l*w_0+E1*sigma+1/2*w_0^2)*k*exp((w_0+l)*S0)+((M_2*k*sigma*E1+1/2*M_2*w_0^2*k+l^2*M_2*k+1/4*w_1*M_4*w_0)*C+1/2*H*k*sigma*variphi*w_0*(R1-1))*(-l*w_0+E1*sigma)*cosh(l*S0))*sinh(l*(-1+S0)))*(l*w_0+E1*sigma))*(l*C*w_1-sigma*E2*R1*k))*(l*C*w_1+sigma*E2*R1*k))*sinh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))-1/2*(w_0^2+4*E1*sigma+4*l^2)^(1/2)*((1/2*(cosh(l)-sinh(l))*(R3*cosh(l*S0)*l+sinh(l*S0)*variphi*sigma)*(-l*w_0+E1*sigma)*(R1-1)*C^2*(l*C*w_1-sigma*E2*R1*k)*(l*w_0+E1*sigma)*M_3*R2*(-1/2*w_1+l)*exp(S0*(w_1-l)-w_1+2*l)+(l*C*w_1+sigma*E2*R1*k)*(1/2*(1/2*w_1+l)*(cosh(l)-sinh(l))*(R3*cosh(l*S0)*l+sinh(l*S0)*variphi*sigma)*(-l*w_0+E1*sigma)*(R1-1)*C^2*(l*w_0+E1*sigma)*M_3*R2*exp(S0*(w_1+l)-w_1)+(1/2*C*k*M_1*(-w_0+l)*(R1-1)*(l*R3*R2*cosh(l*(-1+S0))*R1+variphi*sinh(l*(-1+S0))*sigma)*(-l*w_0+E1*sigma)*exp(-(-w_0+l)*S0)+((-1/2*C*k*R3*M_1*l*(w_0+l)*(R1-1)*exp((w_0+l)*S0)+(-l*w_0+E1*sigma)*(-H*k*R3*l*(R1-1)*cosh(l*S0)+C*(M_2*k*w_0+1/2*M_4*w_1)*sinh(l*S0)))*R1*R2*cosh(l*(-1+S0))-sinh(l*(-1+S0))*(1/2*C*k*M_1*sigma*variphi*(w_0+l)*(R1-1)*exp((w_0+l)*S0)+(-l*w_0+E1*sigma)*cosh(l*S0)*((M_2*k*w_0+1/2*M_4*w_1)*C+H*k*sigma*variphi*(R1-1))))*(l*w_0+E1*sigma))*(l*C*w_1-sigma*E2*R1*k)))*cosh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))+C*k*M_1*w_0*exp(1/2*S0*w_0)*(R1-1)*(l*C*w_1-sigma*E2*R1*k)*(l*C*w_1+sigma*E2*R1*k)*(l*R3*R2*cosh(l*(-1+S0))*R1+variphi*sinh(l*(-1+S0))*sigma)*(E1*sigma+l^2)))*sinh(1/2*(-1+S0)*((w_1^2+4*l^2)*C^2+4*C*sigma*E2*R1*k)^(1/2)/C)-1/8*(-cosh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))*(w_0^2+4*E1*sigma+4*l^2)^(1/2)+sinh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))*w_0)*((w_1^2+4*l^2)*C^2+4*C*sigma*E2*R1*k)^(1/2)*(-l*w_0+E1*sigma)*(l*w_0+E1*sigma)*((-C*M_3*R2*(cosh(l)-sinh(l))*(R1-1)*(R3*cosh(l*S0)*l+sinh(l*S0)*variphi*sigma)*(l*C*w_1-sigma*E2*R1*k)*exp(S0*(w_1-l)-w_1+2*l)+(C*M_3*R2*(cosh(l)-sinh(l))*(R1-1)*(R3*cosh(l*S0)*l+sinh(l*S0)*variphi*sigma)*exp(S0*(w_1+l)-w_1)+2*M_4*(R1*R2*cosh(l*(-1+S0))*sinh(l*S0)-cosh(l*S0)*sinh(l*(-1+S0)))*(l*C*w_1-sigma*E2*R1*k))*(l*C*w_1+sigma*E2*R1*k))*cosh(1/2*(-1+S0)*((w_1^2+4*l^2)*C^2+4*C*sigma*E2*R1*k)^(1/2)/C)-2*C*k*E2*M_3*R1*sigma*R2*exp(1/2*w_1*(-1+S0)+l)*(cosh(l)-sinh(l))*(R1-1)*(R3*cosh(l*S0)*l+sinh(l*S0)*variphi*sigma)))*Phi/H/k/(l^2*C^2*w_1^2-sigma^2*E2^2*R1^2*k^2)/sinh(1/2*(-1+S0)*((w_1^2+4*l^2)*C^2+4*C*sigma*E2*R1*k)^(1/2)/C)/(R1*R2*cosh(l*(-1+S0))*sinh(l*S0)-cosh(l*S0)*sinh(l*(-1+S0)))/(-cosh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))*(w_0^2+4*E1*sigma+4*l^2)^(1/2)+sinh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))*w_0)/(-l^2*w_0^2+E1^2*sigma^2):
 > fsolve({diff((rhs-lhs)(A2),sigma)=0,(rhs-lhs)(A2)=0},{sigma,l=0.5});
 (13)
 > plots:-animate(plot,[[rhs,lhs](A2),sigma=-350..5,-130..10],l=1.5..3.5,frames=20);
 > plots:-animate(plot,[[rhs,lhs](A2),sigma=-350..5,-130..10],l=0..4.5,frames=20);
 > fsolve( [subs(l=4,A2)],sigma=-50, complex);
 >
 (14)
 > fsolve( [subs(l=3.9,A2)],sigma=-51.34547980+1e-5*I, complex);
 (15)
 >
 > fsolve( [subs(l=3.64,A2)],sigma=-22.44118854+1e-5*I, complex);
 (16)
 > eqn2:=Student:-Calculus1:-Roots(subs(l=3.5,A2),sigma);
 (17)
 >
 >
 >
 >

Thanks man. Thats what I was looking for.

But my friend, there is still a problem. In the interval [3.5, 8], sigma has multi real soultions. But the in interval [0, 3.5],  there is a big jump. Now i am looking to find that,  Is there any complex solutions for sigma for `l=0 to 3.5'.

## Re:...

Thanks

I am trying to use a guess value from 3.7 to 8 and to fit a curve then try to extrapolat the rest of the values from 3.7 to 0 backwards. For this I have tried the following code

with(CurveFitting):

Digits:=40;

variphi:=0.38:R1:=0.0009:R2:=8.75:E1:=1:E2:=1:R4:=177.6:C:=1.96:k:=4:H:=7:Pec_i:=-120:Phi:=1:Q:=1.9:R:=20:

S0:=.3592598898850950896270594340253413222287;

F:=1/(R*S0+1-S0);

w_0:=1/k*Pec_i*F*C;

w_1:=Pec_i*F;

M_1:=evalf(Q*Pec_i/k*(H+1/(1-exp(Pec_i))));

M_2:=M_1*exp(w_0*S0);

M_3:=w_1/(1-exp(w_1*(S0-1)));

M_4:=M_3*exp(w_1*(S0-1));

R3:=-38.16;

A:=variphi*sigma*Phi = -2*(((1/4*(cosh(l)-sinh(l))*(R3*cosh(l*S0)*l+sinh(l*S0)*variphi*sigma)*w_0*(-l*w_0+E1*sigma)*(R1-1)*C^2*(l*C*w_1-sigma*E2*R1*k)*(l*w_0+E1*sigma)*M_3*R2*(-1/2*w_1+l)*exp(S0*(w_1-l)-w_1+2*l)+(1/4*(1/2*w_1+l)*(cosh(l)-sinh(l))*(R3*cosh(l*S0)*l+sinh(l*S0)*variphi*sigma)*w_0*(-l*w_0+E1*sigma)*(R1-1)*C^2*(l*w_0+E1*sigma)*M_3*R2*exp(S0*(w_1+l)-w_1)+(-1/2*(l^2-1/2*l*w_0+E1*sigma+1/2*w_0^2)*(-l*w_0+E1*sigma)*(R1-1)*(l*R3*R2*cosh(l*(-1+S0))*R1+variphi*sinh(l*(-1+S0))*sigma)*C*M_1*k*exp(-(-w_0+l)*S0)+((-1/2*l*(R1-1)*C*R3*M_1*(l^2+1/2*l*w_0+E1*sigma+1/2*w_0^2)*k*exp((w_0+l)*S0)+(-l*w_0+E1*sigma)*(-1/2*H*k*R3*l*w_0*(R1-1)*cosh(l*S0)+(M_2*k*sigma*E1+1/2*M_2*w_0^2*k+l^2*M_2*k+1/4*w_1*M_4*w_0)*C*sinh(l*S0)))*R1*R2*cosh(l*(-1+S0))-(1/2*(R1-1)*sigma*variphi*C*M_1*(l^2+1/2*l*w_0+E1*sigma+1/2*w_0^2)*k*exp((w_0+l)*S0)+((M_2*k*sigma*E1+1/2*M_2*w_0^2*k+l^2*M_2*k+1/4*w_1*M_4*w_0)*C+1/2*H*k*sigma*variphi*w_0*(R1-1))*(-l*w_0+E1*sigma)*cosh(l*S0))*sinh(l*(-1+S0)))*(l*w_0+E1*sigma))*(l*C*w_1-sigma*E2*R1*k))*(l*C*w_1+sigma*E2*R1*k))*sinh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))-1/2*(w_0^2+4*E1*sigma+4*l^2)^(1/2)*((1/2*(cosh(l)-sinh(l))*(R3*cosh(l*S0)*l+sinh(l*S0)*variphi*sigma)*(-l*w_0+E1*sigma)*(R1-1)*C^2*(l*C*w_1-sigma*E2*R1*k)*(l*w_0+E1*sigma)*M_3*R2*(-1/2*w_1+l)*exp(S0*(w_1-l)-w_1+2*l)+(l*C*w_1+sigma*E2*R1*k)*(1/2*(1/2*w_1+l)*(cosh(l)-sinh(l))*(R3*cosh(l*S0)*l+sinh(l*S0)*variphi*sigma)*(-l*w_0+E1*sigma)*(R1-1)*C^2*(l*w_0+E1*sigma)*M_3*R2*exp(S0*(w_1+l)-w_1)+(1/2*C*k*M_1*(-w_0+l)*(R1-1)*(l*R3*R2*cosh(l*(-1+S0))*R1+variphi*sinh(l*(-1+S0))*sigma)*(-l*w_0+E1*sigma)*exp(-(-w_0+l)*S0)+((-1/2*C*k*R3*M_1*l*(w_0+l)*(R1-1)*exp((w_0+l)*S0)+(-l*w_0+E1*sigma)*(-H*k*R3*l*(R1-1)*cosh(l*S0)+C*(M_2*k*w_0+1/2*M_4*w_1)*sinh(l*S0)))*R1*R2*cosh(l*(-1+S0))-sinh(l*(-1+S0))*(1/2*C*k*M_1*sigma*variphi*(w_0+l)*(R1-1)*exp((w_0+l)*S0)+(-l*w_0+E1*sigma)*cosh(l*S0)*((M_2*k*w_0+1/2*M_4*w_1)*C+H*k*sigma*variphi*(R1-1))))*(l*w_0+E1*sigma))*(l*C*w_1-sigma*E2*R1*k)))*cosh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))+C*k*M_1*w_0*exp(1/2*S0*w_0)*(R1-1)*(l*C*w_1-sigma*E2*R1*k)*(l*C*w_1+sigma*E2*R1*k)*(l*R3*R2*cosh(l*(-1+S0))*R1+variphi*sinh(l*(-1+S0))*sigma)*(E1*sigma+l^2)))*sinh(1/2*(-1+S0)*((w_1^2+4*l^2)*C^2+4*C*sigma*E2*R1*k)^(1/2)/C)-1/8*(-cosh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))*(w_0^2+4*E1*sigma+4*l^2)^(1/2)+sinh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))*w_0)*((w_1^2+4*l^2)*C^2+4*C*sigma*E2*R1*k)^(1/2)*(-l*w_0+E1*sigma)*(l*w_0+E1*sigma)*((-C*M_3*R2*(cosh(l)-sinh(l))*(R1-1)*(R3*cosh(l*S0)*l+sinh(l*S0)*variphi*sigma)*(l*C*w_1-sigma*E2*R1*k)*exp(S0*(w_1-l)-w_1+2*l)+(C*M_3*R2*(cosh(l)-sinh(l))*(R1-1)*(R3*cosh(l*S0)*l+sinh(l*S0)*variphi*sigma)*exp(S0*(w_1+l)-w_1)+2*M_4*(R1*R2*cosh(l*(-1+S0))*sinh(l*S0)-cosh(l*S0)*sinh(l*(-1+S0)))*(l*C*w_1-sigma*E2*R1*k))*(l*C*w_1+sigma*E2*R1*k))*cosh(1/2*(-1+S0)*((w_1^2+4*l^2)*C^2+4*C*sigma*E2*R1*k)^(1/2)/C)-2*C*k*E2*M_3*R1*sigma*R2*exp(1/2*w_1*(-1+S0)+l)*(cosh(l)-sinh(l))*(R1-1)*(R3*cosh(l*S0)*l+sinh(l*S0)*variphi*sigma)))*Phi/H/k/(l^2*C^2*w_1^2-sigma^2*E2^2*R1^2*k^2)/sinh(1/2*(-1+S0)*((w_1^2+4*l^2)*C^2+4*C*sigma*E2*R1*k)^(1/2)/C)/(R1*R2*cosh(l*(-1+S0))*sinh(l*S0)-cosh(l*S0)*sinh(l*(-1+S0)))/(-cosh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))*(w_0^2+4*E1*sigma+4*l^2)^(1/2)+sinh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))*w_0)/(-l^2*w_0^2+E1^2*sigma^2):

LS:=Matrix(1..81,3):

ss:=0.1;

for j from 38 by 1 to 81 do:

G:=(j-1)*ss;

LS(j,1):=G;

LS(j,2):=simplify(evalf(fsolve(subs(l=G, A)))):

LS(j,3):=LS(j,2):

od:

for j from 1 by 1 to 37 do;

G:=(j-1)*ss:

LS(j,1):=G:

LS(j,2):=ArrayInterpolation(LS(38..80,1),LS(38..80,2), G,degree=20):

od:

for j from 2 by 1 to 37 do;

LS(j,3):=simplify(evalf(fsolve(subs(l=LS(j,1), A)))):

od:

But its also not working I think there is some problem with my approach.

## Re:...

Thanks

I am trying to use a guess value from 3.7 to 8 and to fit a curve then try to extrapolat the rest of the values from 3.7 to 0 backwards. For this I have tried the following code

with(CurveFitting):

Digits:=40;

variphi:=0.38:R1:=0.0009:R2:=8.75:E1:=1:E2:=1:R4:=177.6:C:=1.96:k:=4:H:=7:Pec_i:=-120:Phi:=1:Q:=1.9:R:=20:

S0:=.3592598898850950896270594340253413222287;

F:=1/(R*S0+1-S0);

w_0:=1/k*Pec_i*F*C;

w_1:=Pec_i*F;

M_1:=evalf(Q*Pec_i/k*(H+1/(1-exp(Pec_i))));

M_2:=M_1*exp(w_0*S0);

M_3:=w_1/(1-exp(w_1*(S0-1)));

M_4:=M_3*exp(w_1*(S0-1));

R3:=-38.16;

A:=variphi*sigma*Phi = -2*(((1/4*(cosh(l)-sinh(l))*(R3*cosh(l*S0)*l+sinh(l*S0)*variphi*sigma)*w_0*(-l*w_0+E1*sigma)*(R1-1)*C^2*(l*C*w_1-sigma*E2*R1*k)*(l*w_0+E1*sigma)*M_3*R2*(-1/2*w_1+l)*exp(S0*(w_1-l)-w_1+2*l)+(1/4*(1/2*w_1+l)*(cosh(l)-sinh(l))*(R3*cosh(l*S0)*l+sinh(l*S0)*variphi*sigma)*w_0*(-l*w_0+E1*sigma)*(R1-1)*C^2*(l*w_0+E1*sigma)*M_3*R2*exp(S0*(w_1+l)-w_1)+(-1/2*(l^2-1/2*l*w_0+E1*sigma+1/2*w_0^2)*(-l*w_0+E1*sigma)*(R1-1)*(l*R3*R2*cosh(l*(-1+S0))*R1+variphi*sinh(l*(-1+S0))*sigma)*C*M_1*k*exp(-(-w_0+l)*S0)+((-1/2*l*(R1-1)*C*R3*M_1*(l^2+1/2*l*w_0+E1*sigma+1/2*w_0^2)*k*exp((w_0+l)*S0)+(-l*w_0+E1*sigma)*(-1/2*H*k*R3*l*w_0*(R1-1)*cosh(l*S0)+(M_2*k*sigma*E1+1/2*M_2*w_0^2*k+l^2*M_2*k+1/4*w_1*M_4*w_0)*C*sinh(l*S0)))*R1*R2*cosh(l*(-1+S0))-(1/2*(R1-1)*sigma*variphi*C*M_1*(l^2+1/2*l*w_0+E1*sigma+1/2*w_0^2)*k*exp((w_0+l)*S0)+((M_2*k*sigma*E1+1/2*M_2*w_0^2*k+l^2*M_2*k+1/4*w_1*M_4*w_0)*C+1/2*H*k*sigma*variphi*w_0*(R1-1))*(-l*w_0+E1*sigma)*cosh(l*S0))*sinh(l*(-1+S0)))*(l*w_0+E1*sigma))*(l*C*w_1-sigma*E2*R1*k))*(l*C*w_1+sigma*E2*R1*k))*sinh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))-1/2*(w_0^2+4*E1*sigma+4*l^2)^(1/2)*((1/2*(cosh(l)-sinh(l))*(R3*cosh(l*S0)*l+sinh(l*S0)*variphi*sigma)*(-l*w_0+E1*sigma)*(R1-1)*C^2*(l*C*w_1-sigma*E2*R1*k)*(l*w_0+E1*sigma)*M_3*R2*(-1/2*w_1+l)*exp(S0*(w_1-l)-w_1+2*l)+(l*C*w_1+sigma*E2*R1*k)*(1/2*(1/2*w_1+l)*(cosh(l)-sinh(l))*(R3*cosh(l*S0)*l+sinh(l*S0)*variphi*sigma)*(-l*w_0+E1*sigma)*(R1-1)*C^2*(l*w_0+E1*sigma)*M_3*R2*exp(S0*(w_1+l)-w_1)+(1/2*C*k*M_1*(-w_0+l)*(R1-1)*(l*R3*R2*cosh(l*(-1+S0))*R1+variphi*sinh(l*(-1+S0))*sigma)*(-l*w_0+E1*sigma)*exp(-(-w_0+l)*S0)+((-1/2*C*k*R3*M_1*l*(w_0+l)*(R1-1)*exp((w_0+l)*S0)+(-l*w_0+E1*sigma)*(-H*k*R3*l*(R1-1)*cosh(l*S0)+C*(M_2*k*w_0+1/2*M_4*w_1)*sinh(l*S0)))*R1*R2*cosh(l*(-1+S0))-sinh(l*(-1+S0))*(1/2*C*k*M_1*sigma*variphi*(w_0+l)*(R1-1)*exp((w_0+l)*S0)+(-l*w_0+E1*sigma)*cosh(l*S0)*((M_2*k*w_0+1/2*M_4*w_1)*C+H*k*sigma*variphi*(R1-1))))*(l*w_0+E1*sigma))*(l*C*w_1-sigma*E2*R1*k)))*cosh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))+C*k*M_1*w_0*exp(1/2*S0*w_0)*(R1-1)*(l*C*w_1-sigma*E2*R1*k)*(l*C*w_1+sigma*E2*R1*k)*(l*R3*R2*cosh(l*(-1+S0))*R1+variphi*sinh(l*(-1+S0))*sigma)*(E1*sigma+l^2)))*sinh(1/2*(-1+S0)*((w_1^2+4*l^2)*C^2+4*C*sigma*E2*R1*k)^(1/2)/C)-1/8*(-cosh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))*(w_0^2+4*E1*sigma+4*l^2)^(1/2)+sinh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))*w_0)*((w_1^2+4*l^2)*C^2+4*C*sigma*E2*R1*k)^(1/2)*(-l*w_0+E1*sigma)*(l*w_0+E1*sigma)*((-C*M_3*R2*(cosh(l)-sinh(l))*(R1-1)*(R3*cosh(l*S0)*l+sinh(l*S0)*variphi*sigma)*(l*C*w_1-sigma*E2*R1*k)*exp(S0*(w_1-l)-w_1+2*l)+(C*M_3*R2*(cosh(l)-sinh(l))*(R1-1)*(R3*cosh(l*S0)*l+sinh(l*S0)*variphi*sigma)*exp(S0*(w_1+l)-w_1)+2*M_4*(R1*R2*cosh(l*(-1+S0))*sinh(l*S0)-cosh(l*S0)*sinh(l*(-1+S0)))*(l*C*w_1-sigma*E2*R1*k))*(l*C*w_1+sigma*E2*R1*k))*cosh(1/2*(-1+S0)*((w_1^2+4*l^2)*C^2+4*C*sigma*E2*R1*k)^(1/2)/C)-2*C*k*E2*M_3*R1*sigma*R2*exp(1/2*w_1*(-1+S0)+l)*(cosh(l)-sinh(l))*(R1-1)*(R3*cosh(l*S0)*l+sinh(l*S0)*variphi*sigma)))*Phi/H/k/(l^2*C^2*w_1^2-sigma^2*E2^2*R1^2*k^2)/sinh(1/2*(-1+S0)*((w_1^2+4*l^2)*C^2+4*C*sigma*E2*R1*k)^(1/2)/C)/(R1*R2*cosh(l*(-1+S0))*sinh(l*S0)-cosh(l*S0)*sinh(l*(-1+S0)))/(-cosh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))*(w_0^2+4*E1*sigma+4*l^2)^(1/2)+sinh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))*w_0)/(-l^2*w_0^2+E1^2*sigma^2):

LS:=Matrix(1..81,3):

ss:=0.1;

for j from 38 by 1 to 81 do:

G:=(j-1)*ss;

LS(j,1):=G;

LS(j,2):=simplify(evalf(fsolve(subs(l=G, A)))):

LS(j,3):=LS(j,2):

od:

for j from 1 by 1 to 37 do;

G:=(j-1)*ss:

LS(j,1):=G:

LS(j,2):=ArrayInterpolation(LS(38..80,1),LS(38..80,2), G,degree=20):

od:

for j from 2 by 1 to 37 do;

LS(j,3):=simplify(evalf(fsolve(subs(l=LS(j,1), A)))):

od:

But its also not working I think there is some problem with my approach.

## Re:...

Can we use the value of 3.7 and to get 3.6 untill we reach 0. If we can then how?

3.7 -22.1443585886209766337456856767833538390849200618752824748474
3.8 -19.5885636190466245494729733969928232976910598207760697943391
3.9 -17.5215039974423252872490910655316813866453328121252574492682
4.0 -15.7562118451073083998716254712674127860618363801556087676593
4.1 -14.2035778827336766997300702720157788733317127434423705583457
4.2 -12.8130670646422976067240143018238615607369525904082764850043
4.3 -11.5528805065091562171296493103800204310947315541280249487619
4.4 -10.4015678809169508268921990029382322947525872465781298728833
4.5 -9.34392144075143372387191489011428597848122044232910387364765
4.6 -8.36875018345875378734424554413058044034955522673346676364346
4.7 -7.46757840942950831516874008801141233668471695973444671351081
4.8 -6.63383885626406529496505227845873931631568760262607869628179
4.9 -5.86234847330043717728632579690606776664116618779787344547534
5.0 -5.14895451089142234594739374562417108946677760302629041649626
5.1 -4.49028781633009710194821296840860496356569131552798558204229
5.2 -3.88358612903605971298359068391203573862846688764800130465564
5.3 -3.32656452977122035449005934914122012160013200634417098047139
5.4 -2.81731852219058845710892740702866857729140952599518488519393
5.5 -2.35425023489786118196203489320448236693912219561578313901376
5.6 -1.93601134810602330905403769191000511326251514928889827192753
5.7 -1.56145834346173249235787430975757031137048044916274186181894
5.8 -1.22961698503470400981804828168829172626251435057603745861816
5.9 -.939653818831371541456604399733579630691768068563225354885942
6.0 -.690853080820909477248579854839281120497659887753663449865617
6.1 -.482597824072099173733521725959253752393175776552213060119490
6.2 -.314354374084338817808482770391015053822687232492182609625931
6.3 -.185659436448920637487348387709105932407043557185806108877375
6.4 -.961093380860189427380804591733855189413843503545243162865233e-1
6.5 -.453509995591359605815218931107489985616920865008089972494168e-1
6.6 -.330743230194042948044110457899721797453159321018707301931109e-1
6.7 -.590057462346284979377710812889294830299789880443713008523884e-1
6.8 -.122902763563865218830887491048222028958026732287762061900918
6.9 -.224549253662845888370217531170255221990029403241134850352870
7.0 -.363751484035314895067872409402009704810547115106168118108862
7.1 -.540334686378885150314853771633204210144960186612787216681696
7.2 -.754140115552707501656820661984015240396341802116195497241885
7.3 -1.00502252006136297730989910517873412391303037029861268750648
7.4 -1.29284796406063182039898767228768109701932724962766042989247
7.5 -1.61749195069219310594567656134918206077080821923582046570107
7.6 -1.97883780453755258221328460607676628669795210095136868794865
7.7 -2.37677527752369605228767021140740303399048042216567346988923
7.8 -2.81119934800547161384948812023474349372281734713269216223466
7.9 -3.28200918721917535755560790433157349997360332262564697931579
8.0 -3.78910727102575268195891616911118236824790432001968260476604

## Re:...

Can we use the value of 3.7 and to get 3.6 untill we reach 0. If we can then how?

3.7 -22.1443585886209766337456856767833538390849200618752824748474
3.8 -19.5885636190466245494729733969928232976910598207760697943391
3.9 -17.5215039974423252872490910655316813866453328121252574492682
4.0 -15.7562118451073083998716254712674127860618363801556087676593
4.1 -14.2035778827336766997300702720157788733317127434423705583457
4.2 -12.8130670646422976067240143018238615607369525904082764850043
4.3 -11.5528805065091562171296493103800204310947315541280249487619
4.4 -10.4015678809169508268921990029382322947525872465781298728833
4.5 -9.34392144075143372387191489011428597848122044232910387364765
4.6 -8.36875018345875378734424554413058044034955522673346676364346
4.7 -7.46757840942950831516874008801141233668471695973444671351081
4.8 -6.63383885626406529496505227845873931631568760262607869628179
4.9 -5.86234847330043717728632579690606776664116618779787344547534
5.0 -5.14895451089142234594739374562417108946677760302629041649626
5.1 -4.49028781633009710194821296840860496356569131552798558204229
5.2 -3.88358612903605971298359068391203573862846688764800130465564
5.3 -3.32656452977122035449005934914122012160013200634417098047139
5.4 -2.81731852219058845710892740702866857729140952599518488519393
5.5 -2.35425023489786118196203489320448236693912219561578313901376
5.6 -1.93601134810602330905403769191000511326251514928889827192753
5.7 -1.56145834346173249235787430975757031137048044916274186181894
5.8 -1.22961698503470400981804828168829172626251435057603745861816
5.9 -.939653818831371541456604399733579630691768068563225354885942
6.0 -.690853080820909477248579854839281120497659887753663449865617
6.1 -.482597824072099173733521725959253752393175776552213060119490
6.2 -.314354374084338817808482770391015053822687232492182609625931
6.3 -.185659436448920637487348387709105932407043557185806108877375
6.4 -.961093380860189427380804591733855189413843503545243162865233e-1
6.5 -.453509995591359605815218931107489985616920865008089972494168e-1
6.6 -.330743230194042948044110457899721797453159321018707301931109e-1
6.7 -.590057462346284979377710812889294830299789880443713008523884e-1
6.8 -.122902763563865218830887491048222028958026732287762061900918
6.9 -.224549253662845888370217531170255221990029403241134850352870
7.0 -.363751484035314895067872409402009704810547115106168118108862
7.1 -.540334686378885150314853771633204210144960186612787216681696
7.2 -.754140115552707501656820661984015240396341802116195497241885
7.3 -1.00502252006136297730989910517873412391303037029861268750648
7.4 -1.29284796406063182039898767228768109701932724962766042989247
7.5 -1.61749195069219310594567656134918206077080821923582046570107
7.6 -1.97883780453755258221328460607676628669795210095136868794865
7.7 -2.37677527752369605228767021140740303399048042216567346988923
7.8 -2.81119934800547161384948812023474349372281734713269216223466
7.9 -3.28200918721917535755560790433157349997360332262564697931579
8.0 -3.78910727102575268195891616911118236824790432001968260476604

## re...

 > restart;
 > with(plots):
 > # Based on ZHK's worksheet sent by email 20110202.
 > # Aim is to understand the general shape of the stability problem, especially for very low and high wavenumbers.
 >
 > # Parameter values
 > parvals:= {variphi=0.38, R1=0.0009, R2=8.75, E1=1, E2=1, R4=177.6, C=1.96, k=4, H=7, Pec_i=-120, Phi=1, Q=1.9, R=20};
 (1)
 >
 > # Identifying the possible steady positions of the front
 > frontlocationeq:= (Q*(R*S+1-S)*(H+1/(1-exp(Pec_i)))*exp(C*Pec_i*S/(k*(R*S+1-S)))-H)*(1-exp(Pec_i*(1-S)/(R*S+1-S)))=1;
 (2)
 > Svalues:= Student:-Calculus1:-Roots(subs(parvals,frontlocationeq),S);
 (3)
 > S01:= op(1,Svalues); S02:= op(2,Svalues); S03:= op(3,Svalues);
 (4)
 >
 > # Setting up the dispersion equation
 > F:=1/(R*S0+1-S0):
 > w_0:=1/k*Pec_i*F*C:
 > w_1:=Pec_i*F:
 > M_1:=evalf(Q*Pec_i/k*(H+1/(1-exp(Pec_i)))):
 > M_2:=M_1*exp(w_0*S0):
 > M_3:=w_1/(1-exp(w_1*(S0-1))):
 > M_4:=M_3*exp(w_1*(S0-1)):
 > # The dispersion equation is now A = 0, where
 > A:= -variphi*sigma-2*(((1/4*(R1-1)*(l*w_0+E1*sigma)*R2*(cosh(l)-sinh(l))*(sinh(S0*l)*variphi*sigma+l*R3*cosh(S0*l))*w_0*M_3*(l*C*w_1-sigma*E2*R1*k)*C^2*(-l*w_0+E1*sigma)*(-1/2*w_1+l)*exp(w_1*(S0-1)-l*(S0-2))+(sigma*E2*R1*k+l*C*w_1)*(1/4*(R1-1)*(l*w_0+E1*sigma)*R2*(cosh(l)-sinh(l))*(sinh(S0*l)*variphi*sigma+l*R3*cosh(S0*l))*w_0*M_3*C^2*(1/2*w_1+l)*(-l*w_0+E1*sigma)*exp(w_1*(S0-1)+S0*l)+(l*C*w_1-sigma*E2*R1*k)*(-1/2*(R1-1)*(l^2-1/2*l*w_0+E1*sigma+1/2*w_0^2)*k*M_1*(cosh(l*(S0-1))*R2*R3*l*R1+variphi*sinh(l*(S0-1))*sigma)*C*(-l*w_0+E1*sigma)*exp(-S0*(-w_0+l))+(l*w_0+E1*sigma)*(R2*R1*(-1/2*(R1-1)*l*R3*k*M_1*(l^2+1/2*l*w_0+E1*sigma+1/2*w_0^2)*C*exp(S0*(w_0+l))+(-1/2*H*k*l*w_0*R3*(R1-1)*cosh(S0*l)+C*(k*M_2*sigma*E1+1/2*w_0^2*M_2*k+1/4*w_1*w_0*M_4+l^2*k*M_2)*sinh(S0*l))*(-l*w_0+E1*sigma))*cosh(l*(S0-1))-sinh(l*(S0-1))*(1/2*(R1-1)*k*M_1*(l^2+1/2*l*w_0+E1*sigma+1/2*w_0^2)*C*sigma*variphi*exp(S0*(w_0+l))+cosh(S0*l)*((k*M_2*sigma*E1+1/2*w_0^2*M_2*k+1/4*w_1*w_0*M_4+l^2*k*M_2)*C+1/2*H*k*sigma*w_0*variphi*(R1-1))*(-l*w_0+E1*sigma))))))*sinh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))-1/2*(w_0^2+4*E1*sigma+4*l^2)^(1/2)*((1/2*(R1-1)*(l*w_0+E1*sigma)*R2*(cosh(l)-sinh(l))*(sinh(S0*l)*variphi*sigma+l*R3*cosh(S0*l))*M_3*(l*C*w_1-sigma*E2*R1*k)*C^2*(-l*w_0+E1*sigma)*(-1/2*w_1+l)*exp(w_1*(S0-1)-l*(S0-2))+(sigma*E2*R1*k+l*C*w_1)*(1/2*(R1-1)*(l*w_0+E1*sigma)*R2*(cosh(l)-sinh(l))*(sinh(S0*l)*variphi*sigma+l*R3*cosh(S0*l))*M_3*C^2*(1/2*w_1+l)*(-l*w_0+E1*sigma)*exp(w_1*(S0-1)+S0*l)+(l*C*w_1-sigma*E2*R1*k)*(1/2*C*k*M_1*(-w_0+l)*(R1-1)*(cosh(l*(S0-1))*R2*R3*l*R1+variphi*sinh(l*(S0-1))*sigma)*(-l*w_0+E1*sigma)*exp(-S0*(-w_0+l))+(l*w_0+E1*sigma)*(R2*R1*(-1/2*C*k*l*R3*M_1*(w_0+l)*(R1-1)*exp(S0*(w_0+l))+(-l*w_0+E1*sigma)*(-H*k*l*cosh(S0*l)*R3*(R1-1)+C*(M_2*w_0*k+1/2*M_4*w_1)*sinh(S0*l)))*cosh(l*(S0-1))-sinh(l*(S0-1))*(1/2*C*k*sigma*M_1*variphi*(w_0+l)*(R1-1)*exp(S0*(w_0+l))+cosh(S0*l)*((M_2*w_0*k+1/2*M_4*w_1)*C+H*k*sigma*variphi*(R1-1))*(-l*w_0+E1*sigma))))))*cosh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))+C*k*M_1*w_0*exp(1/2*S0*w_0)*(R1-1)*(l*C*w_1-sigma*E2*R1*k)*(sigma*E2*R1*k+l*C*w_1)*(cosh(l*(S0-1))*R2*R3*l*R1+variphi*sinh(l*(S0-1))*sigma)*(E1*sigma+l^2)))*sinh(1/2*(S0-1)*((w_1^2+4*l^2)*C^2+4*C*sigma*E2*R1*k)^(1/2)/C)-1/8*(l*w_0+E1*sigma)*(-l*w_0+E1*sigma)*((w_1^2+4*l^2)*C^2+4*C*sigma*E2*R1*k)^(1/2)*(-cosh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))*(w_0^2+4*E1*sigma+4*l^2)^(1/2)+sinh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))*w_0)*((-C*R2*M_3*(cosh(l)-sinh(l))*(sinh(S0*l)*variphi*sigma+l*R3*cosh(S0*l))*(R1-1)*(l*C*w_1-sigma*E2*R1*k)*exp(w_1*(S0-1)-l*(S0-2))+(sigma*E2*R1*k+l*C*w_1)*(C*R2*M_3*(cosh(l)-sinh(l))*(sinh(S0*l)*variphi*sigma+l*R3*cosh(S0*l))*(R1-1)*exp(w_1*(S0-1)+S0*l)+2*M_4*(R1*R2*cosh(l*(S0-1))*sinh(S0*l)-cosh(S0*l)*sinh(l*(S0-1)))*(l*C*w_1-sigma*E2*R1*k)))*cosh(1/2*(S0-1)*((w_1^2+4*l^2)*C^2+4*C*sigma*E2*R1*k)^(1/2)/C)-2*C*k*E2*R1*R2*M_3*sigma*exp(1/2*w_1*(S0-1)+l)*(cosh(l)-sinh(l))*(sinh(S0*l)*variphi*sigma+l*R3*cosh(S0*l))*(R1-1)))/H/k/(l^2*C^2*w_1^2-sigma^2*E2^2*R1^2*k^2)/(-w_0^2*l^2+E1^2*sigma^2)/sinh(1/2*(S0-1)*((w_1^2+4*l^2)*C^2+4*C*sigma*E2*R1*k)^(1/2)/C)/(R1*R2*cosh(l*(S0-1))*sinh(S0*l)-cosh(S0*l)*sinh(l*(S0-1)))/(-cosh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))*(w_0^2+4*E1*sigma+4*l^2)^(1/2)+sinh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))*w_0):
 >
 > # To make progress with marginal stability, assume that the exchange of stabilities holds. The condition for marginal stability is now A0 = 0, where
 > A0:= subs(sigma=0,A):
 >
 > # Plot up some cases and see what they look like.
 > Digits:= 40: # High value of Digits required to avoid numerical noise
 > # A0 as a function of R3 and l, for the lowest front position
 > plot3d(subs(S0=S01,subs(parvals,A0)),l=0..10,R3=-500..0,grid=[50,50],axes=BOXED,style=PATCHCONTOUR);
 > display([contourplot(subs(S0=S01,subs(parvals,A0)),l=0..40,R3=-500..0,grid=[50,50],color=black,contours=[0]), contourplot(subs(S0=S01,subs(parvals,A0)),l=0..40,R3=-500..0,grid=[50,50],color=blue)]);
 > p1:=plots:-implicitplot(subs(S0=S01,subs(parvals,A0))=0,l=0..40,R3=-500..0,grid=[50,50]):plots:-display(p1);
 > Q1:=[op([1,1..-2], p1)]:                                                                                                       for i from 1 to nops(Q1) do                                                               fprintf("c:/temp/implotdata.txt","%f",Q1[i]);                                                     fprintf("c:/temp/implotdata.txt","\n");   end do:                                                 fclose("c:/temp/implotdata.txt");
 Error, (in fprintf) number expected for floating point format
 >

## re...

 > restart;
 > with(plots):
 > # Based on ZHK's worksheet sent by email 20110202.
 > # Aim is to understand the general shape of the stability problem, especially for very low and high wavenumbers.
 >
 > # Parameter values
 > parvals:= {variphi=0.38, R1=0.0009, R2=8.75, E1=1, E2=1, R4=177.6, C=1.96, k=4, H=7, Pec_i=-120, Phi=1, Q=1.9, R=20};
 (1)
 >
 > # Identifying the possible steady positions of the front
 > frontlocationeq:= (Q*(R*S+1-S)*(H+1/(1-exp(Pec_i)))*exp(C*Pec_i*S/(k*(R*S+1-S)))-H)*(1-exp(Pec_i*(1-S)/(R*S+1-S)))=1;
 (2)
 > Svalues:= Student:-Calculus1:-Roots(subs(parvals,frontlocationeq),S);
 (3)
 > S01:= op(1,Svalues); S02:= op(2,Svalues); S03:= op(3,Svalues);
 (4)
 >
 > # Setting up the dispersion equation
 > F:=1/(R*S0+1-S0):
 > w_0:=1/k*Pec_i*F*C:
 > w_1:=Pec_i*F:
 > M_1:=evalf(Q*Pec_i/k*(H+1/(1-exp(Pec_i)))):
 > M_2:=M_1*exp(w_0*S0):
 > M_3:=w_1/(1-exp(w_1*(S0-1))):
 > M_4:=M_3*exp(w_1*(S0-1)):
 > # The dispersion equation is now A = 0, where
 > A:= -variphi*sigma-2*(((1/4*(R1-1)*(l*w_0+E1*sigma)*R2*(cosh(l)-sinh(l))*(sinh(S0*l)*variphi*sigma+l*R3*cosh(S0*l))*w_0*M_3*(l*C*w_1-sigma*E2*R1*k)*C^2*(-l*w_0+E1*sigma)*(-1/2*w_1+l)*exp(w_1*(S0-1)-l*(S0-2))+(sigma*E2*R1*k+l*C*w_1)*(1/4*(R1-1)*(l*w_0+E1*sigma)*R2*(cosh(l)-sinh(l))*(sinh(S0*l)*variphi*sigma+l*R3*cosh(S0*l))*w_0*M_3*C^2*(1/2*w_1+l)*(-l*w_0+E1*sigma)*exp(w_1*(S0-1)+S0*l)+(l*C*w_1-sigma*E2*R1*k)*(-1/2*(R1-1)*(l^2-1/2*l*w_0+E1*sigma+1/2*w_0^2)*k*M_1*(cosh(l*(S0-1))*R2*R3*l*R1+variphi*sinh(l*(S0-1))*sigma)*C*(-l*w_0+E1*sigma)*exp(-S0*(-w_0+l))+(l*w_0+E1*sigma)*(R2*R1*(-1/2*(R1-1)*l*R3*k*M_1*(l^2+1/2*l*w_0+E1*sigma+1/2*w_0^2)*C*exp(S0*(w_0+l))+(-1/2*H*k*l*w_0*R3*(R1-1)*cosh(S0*l)+C*(k*M_2*sigma*E1+1/2*w_0^2*M_2*k+1/4*w_1*w_0*M_4+l^2*k*M_2)*sinh(S0*l))*(-l*w_0+E1*sigma))*cosh(l*(S0-1))-sinh(l*(S0-1))*(1/2*(R1-1)*k*M_1*(l^2+1/2*l*w_0+E1*sigma+1/2*w_0^2)*C*sigma*variphi*exp(S0*(w_0+l))+cosh(S0*l)*((k*M_2*sigma*E1+1/2*w_0^2*M_2*k+1/4*w_1*w_0*M_4+l^2*k*M_2)*C+1/2*H*k*sigma*w_0*variphi*(R1-1))*(-l*w_0+E1*sigma))))))*sinh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))-1/2*(w_0^2+4*E1*sigma+4*l^2)^(1/2)*((1/2*(R1-1)*(l*w_0+E1*sigma)*R2*(cosh(l)-sinh(l))*(sinh(S0*l)*variphi*sigma+l*R3*cosh(S0*l))*M_3*(l*C*w_1-sigma*E2*R1*k)*C^2*(-l*w_0+E1*sigma)*(-1/2*w_1+l)*exp(w_1*(S0-1)-l*(S0-2))+(sigma*E2*R1*k+l*C*w_1)*(1/2*(R1-1)*(l*w_0+E1*sigma)*R2*(cosh(l)-sinh(l))*(sinh(S0*l)*variphi*sigma+l*R3*cosh(S0*l))*M_3*C^2*(1/2*w_1+l)*(-l*w_0+E1*sigma)*exp(w_1*(S0-1)+S0*l)+(l*C*w_1-sigma*E2*R1*k)*(1/2*C*k*M_1*(-w_0+l)*(R1-1)*(cosh(l*(S0-1))*R2*R3*l*R1+variphi*sinh(l*(S0-1))*sigma)*(-l*w_0+E1*sigma)*exp(-S0*(-w_0+l))+(l*w_0+E1*sigma)*(R2*R1*(-1/2*C*k*l*R3*M_1*(w_0+l)*(R1-1)*exp(S0*(w_0+l))+(-l*w_0+E1*sigma)*(-H*k*l*cosh(S0*l)*R3*(R1-1)+C*(M_2*w_0*k+1/2*M_4*w_1)*sinh(S0*l)))*cosh(l*(S0-1))-sinh(l*(S0-1))*(1/2*C*k*sigma*M_1*variphi*(w_0+l)*(R1-1)*exp(S0*(w_0+l))+cosh(S0*l)*((M_2*w_0*k+1/2*M_4*w_1)*C+H*k*sigma*variphi*(R1-1))*(-l*w_0+E1*sigma))))))*cosh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))+C*k*M_1*w_0*exp(1/2*S0*w_0)*(R1-1)*(l*C*w_1-sigma*E2*R1*k)*(sigma*E2*R1*k+l*C*w_1)*(cosh(l*(S0-1))*R2*R3*l*R1+variphi*sinh(l*(S0-1))*sigma)*(E1*sigma+l^2)))*sinh(1/2*(S0-1)*((w_1^2+4*l^2)*C^2+4*C*sigma*E2*R1*k)^(1/2)/C)-1/8*(l*w_0+E1*sigma)*(-l*w_0+E1*sigma)*((w_1^2+4*l^2)*C^2+4*C*sigma*E2*R1*k)^(1/2)*(-cosh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))*(w_0^2+4*E1*sigma+4*l^2)^(1/2)+sinh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))*w_0)*((-C*R2*M_3*(cosh(l)-sinh(l))*(sinh(S0*l)*variphi*sigma+l*R3*cosh(S0*l))*(R1-1)*(l*C*w_1-sigma*E2*R1*k)*exp(w_1*(S0-1)-l*(S0-2))+(sigma*E2*R1*k+l*C*w_1)*(C*R2*M_3*(cosh(l)-sinh(l))*(sinh(S0*l)*variphi*sigma+l*R3*cosh(S0*l))*(R1-1)*exp(w_1*(S0-1)+S0*l)+2*M_4*(R1*R2*cosh(l*(S0-1))*sinh(S0*l)-cosh(S0*l)*sinh(l*(S0-1)))*(l*C*w_1-sigma*E2*R1*k)))*cosh(1/2*(S0-1)*((w_1^2+4*l^2)*C^2+4*C*sigma*E2*R1*k)^(1/2)/C)-2*C*k*E2*R1*R2*M_3*sigma*exp(1/2*w_1*(S0-1)+l)*(cosh(l)-sinh(l))*(sinh(S0*l)*variphi*sigma+l*R3*cosh(S0*l))*(R1-1)))/H/k/(l^2*C^2*w_1^2-sigma^2*E2^2*R1^2*k^2)/(-w_0^2*l^2+E1^2*sigma^2)/sinh(1/2*(S0-1)*((w_1^2+4*l^2)*C^2+4*C*sigma*E2*R1*k)^(1/2)/C)/(R1*R2*cosh(l*(S0-1))*sinh(S0*l)-cosh(S0*l)*sinh(l*(S0-1)))/(-cosh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))*(w_0^2+4*E1*sigma+4*l^2)^(1/2)+sinh(1/2*S0*(w_0^2+4*E1*sigma+4*l^2)^(1/2))*w_0):
 >
 > # To make progress with marginal stability, assume that the exchange of stabilities holds. The condition for marginal stability is now A0 = 0, where
 > A0:= subs(sigma=0,A):
 >
 > # Plot up some cases and see what they look like.
 > Digits:= 40: # High value of Digits required to avoid numerical noise
 > # A0 as a function of R3 and l, for the lowest front position
 > plot3d(subs(S0=S01,subs(parvals,A0)),l=0..10,R3=-500..0,grid=[50,50],axes=BOXED,style=PATCHCONTOUR);
 > display([contourplot(subs(S0=S01,subs(parvals,A0)),l=0..40,R3=-500..0,grid=[50,50],color=black,contours=[0]), contourplot(subs(S0=S01,subs(parvals,A0)),l=0..40,R3=-500..0,grid=[50,50],color=blue)]);
 > p1:=plots:-implicitplot(subs(S0=S01,subs(parvals,A0))=0,l=0..40,R3=-500..0,grid=[50,50]):plots:-display(p1);
 > Q1:=[op([1,1..-2], p1)]:                                                                                                       for i from 1 to nops(Q1) do                                                               fprintf("c:/temp/implotdata.txt","%f",Q1[i]);                                                     fprintf("c:/temp/implotdata.txt","\n");   end do:                                                 fclose("c:/temp/implotdata.txt");
 Error, (in fprintf) number expected for floating point format
 >