SetOptions[LogLogPlot,BaseStyle->{FontFamily->"Helvetica",12,GrayLevel[0]},AxesStyle->Black];
SetOptions[Plot,BaseStyle->{FontFamily->"Helvetica",12,GrayLevel[0]},AxesStyle->Black];
SetOptions[ParametricPlot,BaseStyle->{FontFamily->"Helvetica",18,GrayLevel[0]},AxesStyle->Black];

R = 8.314*10^7;   (* universal gas constant; erg/mole/K *)
a13 =2.90*10^17;
a14 = 6.81*10^18;
β13=-1;
β14 = -1;
E13 = 0;                             (* activation energy, erg/mole *)
E14= 496.41*10^10;  (* activation energy, erg/mole *)
T = 5000;                          (* temperature, K *)

r13 = a13*cO[t]*cO[t]*(cO[t]+.35cO2[t])*T^β13*Exp[-E13/R/T];
r14 = a14*cO2[t]*(cO[t]+.35cO2[t])*T^β14*Exp[-E14/R/T];

cOi = 0.001;                 (* initial O concentration mole/cm^3 *)
cO2i = 0.001;              (* initial O2 concentration mole/cm^3 *)
eqns = {
cO'[t] ==-2*r13+2*r14,
cO2'[t]==r13-r14,
cO[0]==cOi,
cO2[0]==cO2i
};
vars = {cO[t],cO2[t]};
tmax = 1*10^(-7);      (* maximum time in seconds *)
sol=NDSolve[eqns,vars,{t,0,tmax}];
ccO =cO[t]/.sol;
ccO2=cO2[t]/.sol;
P =R*T*(ccO+ccO2);  (* ideal gas law for pressure *)
ϵ=10^-11
LogLogPlot[{ccO,ccO2},{t,ϵ,tmax},AxesLabel->{"t (s)","[O],[O2] (mole/cc)"},PlotLegends->{"O","O2"}]
LogLogPlot[P,{t,ϵ,tmax},AxesLabel->{"t (s)","P (dyne/cm2)"},PlotRange->All]
Plot[Abs[ccO+2*ccO2-cOi-2*cO2i],{t,0,tmax},AxesLabel->{"t (s)", "Error (mole/cc)"}]
ccO/.t->tmax
ccO2/.t->tmax

rhs = -2*r13+2*r14/.cO[t]->cccO
f = Expand[rhs/.{cO2[t]->cO2i+cOi/2-cccO/2}]
{sol1,sol2,sol3}=Solve[f==0,cccO]
Plot[f,{cccO,-.0008,.0005},AxesLabel->{"[O] (mole/cc)","f([O])"}]

g=Series[f,{cccO,cccO/.sol3,1}]
τ=-1/Coefficient[g,cccO,1]

sol3
D[f,cccO]/.sol3