Manipulate[ temperature = 25 + 273.15; R = 8.314; Subscript[k, 1][T_] := 1.12 10^3 Exp[-3.12 10^2 (1/T - 1/(25 + 273.15))/R]; Subscript[k, 2][T_] := 2.42 10^3 Exp[-2.34 10^5 (1/T - 1/(25 + 273.15))/R]; Subscript[k, 3][T_] := 6.74 10^2 Exp[-2.64 10^3 (1/T - 1/(25 + 273.15))/R]; Subscript[k, 4][T_] := 1.55 10^3 Exp[-4.94 10^4 (1/T - 1/(25 + 273.15))/R]; Subscript[k, -4][T_] := 1.63 10^(-5) Exp[-1.20 10^2 (1/T - 1/(25 + 273.15))/R]; Subscript[k, 5][T_] := 3.42 10^5 Exp[-1.55 10^3 (1/T - 1/(25 + 273.15))/R]; Subscript[k, 6][T_] := 3.35 10^2 Exp[-9.95 10^2 (1/T - 1/(25 + 273.15))/R]; Subscript[k, a][T_] := 8.09 10^(-4) Exp[-5.73 10^3 (1/T - 1/(25 + 273.15))/R]; Subscript[k, b][T_] := 5.35 10^(-6) Exp[-1.02 10 (1/T - 1/(25 + 273.15))/R]; Subscript[k, app][T_] := 1.13 10^(-1) Exp[-1.39 10^5 (1/T - 1/(25 + 273.15))/R]; Subscript[k, 8][T_] := 1.04 10^(-2) Exp[-1.86 10^3 (1/T - 1/(25 + 273.15))/R]; result0 = NDSolve[{(* governing equations *) D[En[t, temperature], t] == -Subscript[k, 1][temperature] H2O2[t, temperature] En[t, temperature] + Subscript[k, 3][temperature] PhOH[t, temperature] E2[t, temperature] + Subscript[k, a][temperature] E3[t, temperature], D[E1[t, temperature], t] == Subscript[k, 1][temperature] H2O2[t, temperature] En[t, temperature] + (Subscript[k, b][temperature] E3[t, temperature] - Subscript[k, 2][temperature] E1[t, temperature]) PhOH[t, temperature], D[E2[t, temperature], t] == Subscript[k, 2][temperature] PhOH[t, temperature] E1[t, temperature] - (Subscript[k, 3][temperature] PhOH[t, temperature] + Subscript[k, app][temperature] H2O2[t, temperature]) E2[t, temperature], D[E3[t, temperature], t] == Subscript[k, app][temperature] H2O2[t, temperature] E2[t, temperature] - (Subscript[k, a][temperature] + Subscript[k, b][temperature] PhOH[t, temperature]) E3[t, temperature], D[H2O2[t, temperature], t] == -H2O2[t, temperature] (Subscript[k, 1][temperature] En[t, temperature] + Subscript[k, 6][temperature] AmNHOPh[t, temperature] + Subscript[k, app][temperature] E2[t, temperature] + Subscript[k, 8][temperature] H2O2[t, temperature]), D[PhOH[t, temperature], t] == -PhOH[t, temperature] (Subscript[k, 2][temperature] E1[t, temperature] + Subscript[k, 3][temperature] E2[t, temperature] + Subscript[k, -4][temperature] AmNH[t, temperature] + Subscript[k, b][temperature] E3[t, temperature]) + Subscript[k, 4][temperature] PhO[t, temperature] AmNH2[t, temperature], D[AmNH2[t, temperature], t] == -Subscript[k, 4][temperature] AmNH2[t, temperature] PhO[t, temperature] + Subscript[k, -4][temperature] AmNH[t, temperature] PhOH[t, temperature], D[PhO[t, temperature], t] == (Subscript[k, 2][temperature] E1[t, temperature] + Subscript[k, 3][temperature] E2[t, temperature] + Subscript[k, -4][temperature] AmNH[t, temperature] + Subscript[k, b][temperature] E3[t, temperature]) PhOH[t, temperature] - (Subscript[k, 4][temperature] AmNH2[t, temperature] + Subscript[k, 5][temperature] AmNH[t, temperature]) PhO[t, temperature], D[AmNH[t, temperature], t] == Subscript[k, 4][temperature] PhO[t, temperature] AmNH2[t, temperature] - (Subscript[k, -4][temperature] PhOH[t, temperature] + Subscript[k, 5][temperature] PhO[t, temperature]) AmNH[t, temperature], D[AmNHOPh[t, temperature], t] == Subscript[k, 5][temperature] PhO[t, temperature] AmNH[t, temperature] - Subscript[k, 6][temperature] AmNHOPh[t, temperature] H2O2[t, temperature], D[Dye[t, temperature], t] == Subscript[k, 6][temperature] AmNHOPh[t, temperature] H2O2[t, temperature], (* E[t,temperature]==E[0,temperature]-E1[t,temperature]-E2[t, temperature]-E3[t,temperature], E1[t,temperature]==(Subscript[k, 3][temperature]E2[t, temperature]+(Subscript[k, b][temperature]+Subscript[k, a][ temperature]/PhOH[t,temperature])E3[t,temperature])/Subscript[k, 2][temperature], E2[t,temperature]==(E[0,temperature]-((Subscript[k, b][ temperature]PhOH[t,temperature]+Subscript[k, a][ temperature])/(Subscript[k, 2][temperature]PhOH[t,temperature])+ Subscript[k, a][temperature]/(Subscript[k, 1][temperature]H2O2[t, temperature])+1)E3[t,temperature])/(Subscript[k, 3][temperature]/ Subscript[k, 2][temperature]+(Subscript[k, 3][temperature]PhOH[t, temperature])/(Subscript[k, 1][temperature]H2O2[t,temperature])+1), E3[t,temperature]==(E[0,temperature]Subscript[k, app][ temperature])/(((Subscript[k, b][temperature]PhOH[t,temperature]+ Subscript[k, a][temperature])/H2O2[t,temperature])(Subscript[k, 3][temperature]/Subscript[k, 2][temperature]+(Subscript[k, 3][ temperature]PhOH[t,temperature])/(Subscript[k, 1][ temperature]H2O2[t,temperature])+1)+((Subscript[k, b][ temperature]PhOH[t,temperature]+Subscript[k, a][temperature])/( Subscript[k, 2][temperature]PhOH[t,temperature])+Subscript[k, a][ temperature]/(Subscript[k, 1][temperature]H2O2[t, temperature])+1)Subscript[k, app][temperature]), PhO[t,temperature]==(Subscript[k, 2][temperature]E1[t, temperature]+Subscript[k, 3][temperature]E2[t,temperature]+ Subscript[k, -4][temperature]AmNH[t,temperature]+Subscript[k, b][ temperature]E3[t,temperature])PhOH[t,temperature]/(Subscript[k, 4][temperature]AmNH2[t,temperature]+Subscript[k, 5][ temperature]AmNH[t,temperature]), *) (* initial conditions *) En[0, temperature] == 0.1/(40 10^3), E1[0, temperature] == 0, E2[0, temperature] == 0, E3[0, temperature] == 0, H2O2[0, temperature] == 0.1, PhOH[0, temperature] == 25, AmNH2[0, temperature] == 0.4, PhO[0, temperature] == 0, AmNH[0, temperature] == 0, AmNHOPh[0, temperature] == 0, Dye[0, temperature] == 0 }, {En[t, temperature], E1[t, temperature], E2[t, temperature], E3[t, temperature], H2O2[t, temperature], PhOH[t, temperature], AmNH2[t, temperature], PhO[t, temperature], AmNH[t, temperature], AmNHOPh[t, temperature], Dye[t, temperature]}, {t, 0, 600}]; temperature = temperature0 + 273.15; result = NDSolve[{(* governing equations *) D[En[t, temperature], t] == -Subscript[k, 1][temperature] H2O2[t, temperature] En[t, temperature] + Subscript[k, 3][temperature] PhOH[t, temperature] E2[t, temperature] + Subscript[k, a][temperature] E3[t, temperature], D[E1[t, temperature], t] == Subscript[k, 1][temperature] H2O2[t, temperature] En[t, temperature] + (Subscript[k, b][temperature] E3[t, temperature] - Subscript[k, 2][temperature] E1[t, temperature]) PhOH[t, temperature], D[E2[t, temperature], t] == Subscript[k, 2][temperature] PhOH[t, temperature] E1[t, temperature] - (Subscript[k, 3][temperature] PhOH[t, temperature] + Subscript[k, app][temperature] H2O2[t, temperature]) E2[t, temperature], D[E3[t, temperature], t] == Subscript[k, app][temperature] H2O2[t, temperature] E2[t, temperature] - (Subscript[k, a][temperature] + Subscript[k, b][temperature] PhOH[t, temperature]) E3[t, temperature], D[H2O2[t, temperature], t] == -H2O2[t, temperature] (Subscript[k, 1][temperature] En[t, temperature] + Subscript[k, 6][temperature] AmNHOPh[t, temperature] + Subscript[k, app][temperature] E2[t, temperature] + Subscript[k, 8][temperature] H2O2[t, temperature]), D[PhOH[t, temperature], t] == -PhOH[t, temperature] (Subscript[k, 2][temperature] E1[t, temperature] + Subscript[k, 3][temperature] E2[t, temperature] + Subscript[k, -4][temperature] AmNH[t, temperature] + Subscript[k, b][temperature] E3[t, temperature]) + Subscript[k, 4][temperature] PhO[t, temperature] AmNH2[t, temperature], D[AmNH2[t, temperature], t] == -Subscript[k, 4][temperature] AmNH2[t, temperature] PhO[t, temperature] + Subscript[k, -4][temperature] AmNH[t, temperature] PhOH[t, temperature], D[PhO[t, temperature], t] == (Subscript[k, 2][temperature] E1[t, temperature] + Subscript[k, 3][temperature] E2[t, temperature] + Subscript[k, -4][temperature] AmNH[t, temperature] + Subscript[k, b][temperature] E3[t, temperature]) PhOH[t, temperature] - (Subscript[k, 4][temperature] AmNH2[t, temperature] + Subscript[k, 5][temperature] AmNH[t, temperature]) PhO[t, temperature], D[AmNH[t, temperature], t] == Subscript[k, 4][temperature] PhO[t, temperature] AmNH2[t, temperature] - (Subscript[k, -4][temperature] PhOH[t, temperature] + Subscript[k, 5][temperature] PhO[t, temperature]) AmNH[t, temperature], D[AmNHOPh[t, temperature], t] == Subscript[k, 5][temperature] PhO[t, temperature] AmNH[t, temperature] - Subscript[k, 6][temperature] AmNHOPh[t, temperature] H2O2[t, temperature], D[Dye[t, temperature], t] == Subscript[k, 6][temperature] AmNHOPh[t, temperature] H2O2[t, temperature], (* E[t,temperature]==E[0,temperature]-E1[t,temperature]-E2[t, temperature]-E3[t,temperature], E1[t,temperature]==(Subscript[k, 3][temperature]E2[t, temperature]+(Subscript[k, b][temperature]+Subscript[k, a][ temperature]/PhOH[t,temperature])E3[t,temperature])/Subscript[k, 2][temperature], E2[t,temperature]==(E[0,temperature]-((Subscript[k, b][ temperature]PhOH[t,temperature]+Subscript[k, a][ temperature])/(Subscript[k, 2][temperature]PhOH[t,temperature])+ Subscript[k, a][temperature]/(Subscript[k, 1][temperature]H2O2[t, temperature])+1)E3[t,temperature])/(Subscript[k, 3][temperature]/ Subscript[k, 2][temperature]+(Subscript[k, 3][temperature]PhOH[t, temperature])/(Subscript[k, 1][temperature]H2O2[t,temperature])+1), E3[t,temperature]==(E[0,temperature]Subscript[k, app][ temperature])/(((Subscript[k, b][temperature]PhOH[t,temperature]+ Subscript[k, a][temperature])/H2O2[t,temperature])(Subscript[k, 3][temperature]/Subscript[k, 2][temperature]+(Subscript[k, 3][ temperature]PhOH[t,temperature])/(Subscript[k, 1][ temperature]H2O2[t,temperature])+1)+((Subscript[k, b][ temperature]PhOH[t,temperature]+Subscript[k, a][temperature])/( Subscript[k, 2][temperature]PhOH[t,temperature])+Subscript[k, a][ temperature]/(Subscript[k, 1][temperature]H2O2[t, temperature])+1)Subscript[k, app][temperature]), PhO[t,temperature]==(Subscript[k, 2][temperature]E1[t, temperature]+Subscript[k, 3][temperature]E2[t,temperature]+ Subscript[k, -4][temperature]AmNH[t,temperature]+Subscript[k, b][ temperature]E3[t,temperature])PhOH[t,temperature]/(Subscript[k, 4][temperature]AmNH2[t,temperature]+Subscript[k, 5][ temperature]AmNH[t,temperature]), *) (* initial conditions *) En[0, temperature] == concentration/(40 10^3), E1[0, temperature] == 0, E2[0, temperature] == 0, E3[0, temperature] == 0, H2O2[0, temperature] == peroxide, PhOH[0, temperature] == 25, AmNH2[0, temperature] == 0.4, PhO[0, temperature] == 0, AmNH[0, temperature] == 0, AmNHOPh[0, temperature] == 0, Dye[0, temperature] == 0 }, {En[t, temperature], E1[t, temperature], E2[t, temperature], E3[t, temperature], H2O2[t, temperature], PhOH[t, temperature], AmNH2[t, temperature], PhO[t, temperature], AmNH[t, temperature], AmNHOPh[t, temperature], Dye[t, temperature]}, {t, 0, 600}]; Plot[{Dye[t, temperature] /. result, Dye[t, 25 + 273.15] /. result0}, {t, 0, 600}, PlotRange -> Full, PlotStyle -> {{Blue}, {Red}}, AxesLabel -> {"time (s)", "Dye (mM)"}], {{temperature0, 25, "temperature in C"}, 10, 40, 1}, Delimiter, {{concentration, 0.1, "HRP in mg/L"}, 0.03, 0.6}, {{peroxide, 0.1, "H2O2 in mM"}, 0.03, 0.2}]