This Workbook sets up, simulates and graphs the  equations that describe the reaction chain in section 3.4 of SB2:
x1 <=> x2 <=> x3 --> x4
BOP Sept  15, 2007
Addiel de Alba. April,2010

Step 1: define the equations

rhs1=b1-v1;
rhs2=v1-v2;
rhs3=v2-v3;
rhs4=v3;

step 2 : define the parmeter values and rate laws

k1=1;k2=0.01;k3=0.0001;
K1=1.0;K2=1.0;
v1=k1*(x1[t]-x2[t]/K1);
v2=k2*(x2[t]-x3[t]/K2);
v3=k3*x3[t];

step 3:  set up the conditions and simulate

b1=0;
tfinal=50000;
solution=NDSolve[{x1'[t]==rhs1,x2'[t]==rhs2,x3'[t]==rhs3,x4'[t]==rhs4,x1[0]==1.0,x2[0]==0,x3[0]==0,x4[0]==0.},{x1,x2,x3,x4},{t,0,tfinal}];

step 4:  graph the results

<<"PlotLegends`"
graph1=Plot[{x1[t]/.solution,x2[t]/.solution,x3[t]/.solution,x4[t]/.solution},{t,0.01`,3},PlotLabel->"Concentrations vs time",PlotLegend->{"x1","x2","x3","x4"},AxesLabel->{"Time"," "},LegendPosition->{.6`,0.2`},PlotStyle->{Hue[0.5`],Hue[0.6`],Hue[0.8`],Hue[1.`]},PlotRange->{0,1.`},LegendSize->{0.35`,0.35`},ImageSize->{600}]
graph2=Plot[{x1[t]/.solution,x2[t]/.solution,x3[t]/.solution,x4[t]/.solution},{t,3,300},PlotLabel->"Concentrations vs time",PlotLegend->{"x1","x2","x3","x4"},AxesLabel->{"Time"," "},LegendPosition->{.6`,0.2`},PlotStyle->{Hue[0.5`],Hue[0.6`],Hue[0.8`],Hue[1.`]},PlotRange->{0,1.`},LegendSize->{0.35`,0.35`},ImageSize->{600}]
graph3=Plot[{x1[t]/.solution,x2[t]/.solution,x3[t]/.solution,x4[t]/.solution},{t,300,tfinal},PlotLabel->"Concentrations vs time",PlotLegend->{"x1","x2","x3","x4"},AxesLabel->{"Time"," "},LegendPosition->{.6`,0.2`},PlotStyle->{Hue[0.5`],Hue[0.6`],Hue[0.8`],Hue[1.`]},PlotRange->{0,1.`},LegendSize->{0.35`,0.35`},ImageSize->{600}]
graph4=ParametricPlot[Evaluate[{x1[t],x2[t]}/.solution],{t,0,tfinal},PlotLabel->"Dynamic Phase portrait",PlotStyle->{Hue[0.5`]},AxesLabel->{"x1","x2"},PlotRange->{{0,1},{0,0.5`}},ImageSize->{600}]

Step 5:  5xport the graphs

Export["Section3.4A.eps",graph1];
Export["Section3.4B.eps",graph2];
Export["Section3.4C.eps",graph3];
Export["Section3.4D.eps",graph4];