Computing multiple steady states 
in a nonisothermal CSTR
Author's Data: Housam BINOUS
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Acknowledgement :

This example of the was given by Kenneth J. Beers in Numerical Methods for Chemical Engineering, Applications in Matlab, Cambridge University Press, 2007.

At high activation energies (γ=12), we find multiple steady states.

Off[Solve::"ifun"]

χ=0;θc=1;β=-1;γ=12;

Solve[1-θ-β 10^Da Exp[γ (θ-1)/(θ)] φA -χ (θ-θc)==0,φA]

Solve[1-φA -10^Da Exp[γ (θ-1)/(θ)] φA==0,θ]

1-φA -10^Da Exp[γ (θ-1)/(θ)] φA==0 /.φA->(10^(-Da)) (E^(-12+(12/θ))) (-1+θ)

plt1=ContourPlot[1-(10^(-Da)) (E^(-12+(12/θ))) (-1+θ)-(E^(-12+(12/θ)+((12 (-1+θ))/θ))) (-1+θ)==0,
{Da,-2,0},{θ,1,2},ContourStyle->Hue[0.5]]

1-θ-β 10^Da Exp[γ (θ-1)/(θ)] φA -χ (θ-θc)==0/.θ->-(12/(-12+Log[(((10^(-Da)) (1-φA))/(φA))]))

plt2=ContourPlot[1+(10^Da) (E^((-1-(12/(-12+Log[(((10^(-Da)) (1-φA))/(φA))]))) (12-Log[(((10^(-Da)) (1-φA))/(φA))]))) φA+(12/(-12+Log[(((10^(-Da)) (1-φA))/(φA))]))==0,
{Da,-2,0},{φA,0.001,0.999},ContourStyle->Hue[0.8]]

With ImplicitPlot, we get the same figure as figure 4.20, page 207, given in the book by Kenneth J. Beers in Numerical Methods for Chemical Engineering, Applications in Matlab, Cambridge University Press, 2007.

Show[GraphicsArray[{{plt2},{plt1}}]]