Numerical study of an oscillatory chemical reaction

Or: Oscillations beyond the harmonic oscillator

The Notebook courtesy of Markus Deserno.	

Please note:
(1) In order to evaluate a line, you need to press SHIFT+RETURN.
(2) If you change parameters, you need to re-evaluate the subsequent expressions as well or run Evaluate Notebook.

Here are the two parameters for the differential equation:
(Notice that this script, for simplicity, calls them a and b rather than α and β, and that the scaled components are called x and y rather than ξ and η.)

a := 1;

b := 3;

Here's the differential equation, i.e., the vector field of the flow

Vx[x_,y_] := 1 + a x^2 y - (b+1)x

Vy[x_,y_] := -a x^2 y + b x

Here's a plot of that vector field (notice that the length of the vectors is scaled by Log[4+vectorlength], so that we can see them better).

VectorPlot[{Vx[x,y],Vy[x,y]},{x,0,5},{y,0,5},PlotRange->{{0,5},{0,5}},VectorScale->{Automatic,Automatic,Log[4+#5]&}]

Here' s an illustration of the flow-field belonging to that vector field

StreamPlot[{Vx[x,y],Vy[x,y]},{x,0,5},{y,0,5},StreamPoints->Fine,PlotRange->{{0,5},{0,5}}]

Here's the initial condition for the differential equation, i.e., from what values the two (scaled) concentrations x and y start at t=0.  You should play around with where you start to integrate.

x0 := 1.1;

y0 := 3;

This line solves the differential equation

solution=NDSolve[{X'[t]==Vx[X[t],Y[t]],Y'[t]==Vy[X[t],Y[t]],X[0]==x0,Y[0]==y0},{X,Y},{t,0,100}];

Here's a plot of the variable x versus time

Plot[Evaluate[{X[t]}/.solution],{t,0,50}]

Here's a plot of the variable y versus time

Plot[Evaluate[{Y[t]}/.solution],{t,0,50}]

Here's a parametric plot of the trajectory (x(t),y(t)) as a function of time.  Also try starting the integration from different initial points (e.g., close to the stationary point or far away) and see what you get.  The plot also includes the so-called "null-clines", which are the lines where either the vertical component of the vector field vanishes and thus the solution is horizontal (red) or the horizontal component of the vector field vanishes and thus the solution is vertical (blue).  Notice that the two null-clines of course intersect in the stationary point.

ParametricPlot[{Evaluate[{X[t],Y[t]}/.solution],{t,(b/(a t))},{t,(((b+1) t-1)/(a (t^2)))}},{t,0,100},PlotRange->{{0,5},{0,5}},AspectRatio->1,PlotPoints->500,PlotStyle->{{Thickness[0.005]},{Thickness[0.003],RGBColor[1,0,0]},{Thickness[0.003],RGBColor[0,0,1]}}]