restart:with(MultiSeries):with(plots):
parvalues:={C=2.02,H_liq=8.74,R=22.11,R1=0.0006};
G:=R3*k*(1-R1)/(C*Pec_i);
F1:=1/(R*S+1-S);
F:=F1*(1-R*G*S);
F2:=Pec_i*F;
A1:=Theta0/C=(exp(F2*(1-S))-1)*(H_liq+1/(1-exp(-C/k*F2*(S))));
R3values:= [0,-1,-3,-4]; colours:= [red,blue,green,black]; display([seq(plottools[transform]( unapply([k, subs(parvalues, Pec_i = -20, R3=R3values[i], F2)], (k,S))) (implicitplot(subs(parvalues,Theta0=4,R3=R3values[i...