gasConst1 = 1.987;
gasConst2 = 0.08206;
gasConst3 = 8.314;
avogadro = 6.022*10^23;
temp = 973;
pZero = 1.0;
c2h6Zero = pZero/gasConst2/temp
a1 = 4.0*10^16;
a2 = 3.8*10^11;
a3 = 3.2*10^13;
a4 = 1.0*10^11;
a5 = 4.0*10^10;
k1 = 
        a1*Exp[-87500/gasConst1/temp];
k2 = 
        a2*Exp[-16500/gasConst1/temp];
k3 = 
        a3*Exp[-40000/gasConst1/temp];
k4 = 
        a4*Exp[-9700/gasConst1/temp];
k5 = a5*Exp[-0/gasConst1/temp];

concCH3 = 2*k1/k2 // N
concC2H5 = (k1*k4/(k3*k5))^0.5*
        concC2H6
concH = (k1*k3/(k4*k5))^0.5

c2h6 = 
      DSolve[{-y1'[t] == (2*k1 + (k1*k3*k4/k5)^0.5)*
              y1[t], y1[0] == c2h6Zero}, y1[t], t]

Plot[y1[t] /. c2h6, {t, 0, 10}];

c2h4 = 
      DSolve[{y2'[t] == (k1*k3*k4/k5)^0.5*y1[t] /. c2h6, 
          y2[0] == 0}, y2[t], t]

ch4 = 
      DSolve[{y3'[t] == k1*y1[t] /. c2h6, y3[0] == 0}, y3[t], 
        t]

h2 = 
      DSolve[{y4'[t] == (k1*k3*k4/k5)^0.5*y1[t] /. c2h6, 
          y4[0] == 0}, y4[t], t]

Plot[y2[t] /. c2h4, {t, 0, 10}];

Plot[y3[t] /. ch4, {t, 0, 10}];

Plot[y4[t] /. h2, {t, 0, 10}];

Now, of course, comes the most interesting part of the analysis: (a) the effect of temperature and pressure on the product composition; (b) the sensitivity of the results to the kinetic parameters (measured or assumed).

LRR3@psu.edu (modified 10/25/04)