The Semenov Model
Author's Data: Housam Binous
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Details about the Semenov model can be found in: 
G. Nicolis, Dynamique Chimique, Thermodynamique, cinétique et mécanique statistique, Dunod, Paris 2005. 
The author would like to acknowledge the help of Professor
Brian Higgins, Chemical Engineering Department, U C Davis who kindly provided notebooks on continuation methods.

The Semenov model can be written as follows:

		(dT/d)τ = 
        exp(-(1/(ϵ T))) - γ (T - 1)

Let us solve the non linear ODE for ϵ=0.1 and γ=0.5.

Off[General::"spell1"]

ϵ = 0.1; γ = 0.5;

sol = 
      NDSolve[{D[T[τ], {τ, 1}] ==  
            Exp[-1/(ϵ T[τ])] - γ (T[τ] - 1), T[0] ==  10}, {T[τ]}, {τ, 0, 1000}]

Plot[T[τ] /. sol, {τ, 0, 100}, 
      PlotStyle ->  RGBColor[1, 0, 0], 
      PlotRange ->  All]

T[τ] /. sol /. τ ->  1000

Let us find the steady state of the non linear ODE as a function of γ for ϵ=0.1 using ImplicitPlot.

γ =.

Needs["Graphics`ImplicitPlot`"]

plt2 = 
      ImplicitPlot[
        Exp[-1/(ϵ T)] - γ (T - 1) == 
          0 , {γ, 0, 0.5}, {T, 0.01, 100}, AspectRatio ->  .7, 
        PlotStyle -> RGBColor[1, 0, 0], 
        PlotPoints ->  1000]

Let us find the steady state of the non linear ODE as a function of γ for ϵ=0.1 using the continuation method. Indeed, we use the arc length continuation method to solve a single equation f1 (T, γ) = 0 that exhibited  turning points.

f1 = 
        Exp[-1/(ϵ T[s])] - γ[
              s] (T[s] - 1);

Outer[D, {f1}, {γ[s], s}]

Solve[%[[1]].{γ'[s],1}== - First[Outer[D, {f1}, {T[s]}][[1]] T'[s]],γ'[s]]//Simplify

eq1=γ'[s]== ((\(E\^\(-\(10.`\/T[s]\)\)\) \((\(-20.`\) + 
                2.`\ E\^\(10.`\/T[s]\)\ T[s]\^2\ γ[s])\) T'[s])/((2.` - 
                2.` T[s]) T[s]^2))

eq2=  T'[s]== \(\(+1\)/\@\((1 + \((\(\(E\^\(-\(10.`\/T[s]\)\)\)\(\ \)\((\(-20.`\) + 2.`\ E\^\(10.`\/T[s]\)\ T[s]\^2\ γ[s])\)\(\ \)\)\/\(\((\(\(2.`\)\(\)\) - 2.`\ T[s])\)\ T[s]\^2\))\)^2)\)\)

Sol = 
      NDSolve[{eq1, eq2, 
            T[0] == 1.00009`, γ[0] == 0.5`}, {T[
              s], γ[s]}, {s, 0, 100}, MaxSteps -> 1000] // 
        Timing

plt4 = 
        ParametricPlot[
          Evaluate[{γ[s], T[s]} /. Sol[[2]]], {s, 0, 100}, 
          PlotRange -> All, 
          PlotStyle -> RGBColor[0, 0, 1]];

The arc length continuation method and the method based on ImplicitPlot give the same result. We get the well known S-shape.

Show[plt2, plt4]