# Return to Question

 Post Closed as "off-topic" by user9660, m_goldberg, MarcoB, Yves Klett, Öskå occurred Dec 24 '15 at 10:44 4 added 531 characters in body edited May 9 '15 at 15:32 André Oliveira 5355 bronze badges  Tsup[Φ_] := To + ΔT * Cos[Φ]; Rsup[Φ_] := Ro + ΔR * Cos[Φ + π]; To=2200 Ro=2.1 ΔT=440 ΔR=0.25  Having t[r_,Φ_] for each 'r', I can calculate 'k', wich goes from k1 to k143, to each t[r_,Φ_]: ki = α*(t[r,Φ]/300)^βExp[-ɣ * Exp[-ɣ/t[r,Φ]] Then I can finally calcute the molecular concentration profiles by solving a system of 50 non-linear, stiff, ordinary and coupled differential equations. Here an example with three of those equations. nH'[t] ==, ((k17*nH2[t]*nO[t]where +α,β k18*nH2[t]*nOH[t]and +ɣ k21*nC[t]*nOH[t]depend +on  the species. In the example k62*nO[t]*nOH[t]above, +these k97*nOH[t]*nCO[t]are +the  'k' to solve the system k141*nTi[t]*nOH[t])of ODE: k1=(6.99*10^- 14)*(k1*nH[t]*nOH[t] + k2*nH[t]*nH2O[t] + k5*nH[t]*nCO[t] + k8*nH[t]*nO2[t] + k10*nH[t]*nCO2[t](t[r,Φ]/300)^2.8)*Exp[-1950/t[r, Φ]]; nH2'[t] == ((k1*nH[t]*nOH[t] + k2*nH[t]*nH2O[t] +k2=2*10^-10 k143*nTi[t]*nH2O[t]for +all temperatures; k144*nTiO[t]*nH2O[t]) k17=(3.14*10^- 13)*(k17*nH2[t]*nO[t] + k18*nH2[t]*nOH[t] + k19*nH2[t]*nO2[t](t[r,Φ]/300)^2.7)*Exp[-3150/t[r, Φ]]; nC'[t] == (k18=(k5*nH[t]*nCO[t] + k140*nTi[t]*nCO[t]) 2.05*10^- 12)*(k20*nC[t]*nOH[t] + k21*nC[t]*nOH[t] + k22*nC[t]*nCO[t] + k26*nC[t]*nO2[t](t[r,Φ]/300)^1.52)*Exp[-1736/t[r, Φ]]; ..k62=(1.77*10^-11)*Exp[178/t[r, Φ]]; nH[0] = 10^k63=(1.85*10^-711)*((t[r, nH2[0] = 10^Φ]/300)^0.95)*Exp[-28571/t[r, Φ]]; nC[0]k94 = 10^(1.65*10^-312)*((t[r, ..Φ]/300)^1.14)*Exp[-50/t[r,Φ]]  Then I can finally calcute the molecular concentration profiles by solving a system of 50 non-linear, stiff, ordinary and coupled differential equations. Here an example of those equations.  s=NDSolve[{ nH'[t] == ((k17*nH2[t]*nO[t] + k18*nH2[t]*nOH[t] + k62*nO[t]*nOH[t]) - (k1*nH[t]*nOH[t] + k2*nH[t]*nH2O[t])), nH2'[t] == ((k1*nH[t]*nOH[t] + k2*nH[t]*nH2O[t]) - (k17*nH2[t]*nO[t] + k18*nH2[t]*nOH[t])), nO'[t] == ((k1*nH[t]*nOH[t] + k94*nOH[t]*nOH[t]) - (k62*nO[t]*nOH[t] + k63*nO[t]*nH2O[t])), nOH'[t] == ((k2*nH[t]*nH2O[t] + k17*nH2[t]*nO[t] + k63*nO[t]*nH2O[t]) - (k1*nH[t]*nOH[t] + k18*nH2[t]*nOH[t] + k62*nO[t]*nOH[t] + k94*nOH[t]*nOH[t])), nH2O'[t] == ((k18*nH2[t]*nOH[t] + k94*nOH[t]*nOH[t]) - (k2*nH[t]*nH2O[t] + k63*nO[t]*nH2O[t])), ..., nH[0] == 10^-7, nH2[0] == 5*10^-1, nO[0] == 10^-5, nOH[0] == 4.9*10^-7, nH2O[0] == 4*10^-6, }, ] When solving this system, the program will give me the molecular concentration for each specie which I'd use as initial condition to the next iterationinitial condition to the next iteration ; After solving the system, the program should go back to step 4, use the next value of t[r_,Φ_], calculate k1 to k143 again and use them to solve the system using the previously molecular concentration (nX) that was calculated before, until all values of t[r_,Φ_] are used.  Tsup[Φ_] := To + ΔT * Cos[Φ]; Rsup[Φ_] := Ro + ΔR * Cos[Φ + π];   Having t[r_,Φ_] for each 'r', I can calculate 'k', wich goes from k1 to k143, to each t[r_,Φ_]: ki = α*(t[r,Φ]/300)^βExp[-ɣt[r,Φ]] Then I can finally calcute the molecular concentration profiles by solving a system of 50 non-linear, stiff, ordinary and coupled differential equations. Here an example with three of those equations. nH'[t] == ((k17*nH2[t]*nO[t] + k18*nH2[t]*nOH[t] + k21*nC[t]*nOH[t] +   k62*nO[t]*nOH[t] + k97*nOH[t]*nCO[t] +   k141*nTi[t]*nOH[t]) - (k1*nH[t]*nOH[t] + k2*nH[t]*nH2O[t] + k5*nH[t]*nCO[t] + k8*nH[t]*nO2[t] + k10*nH[t]*nCO2[t])), nH2'[t] == ((k1*nH[t]*nOH[t] + k2*nH[t]*nH2O[t] + k143*nTi[t]*nH2O[t] + k144*nTiO[t]*nH2O[t]) - (k17*nH2[t]*nO[t] + k18*nH2[t]*nOH[t] + k19*nH2[t]*nO2[t])), nC'[t] == ((k5*nH[t]*nCO[t] + k140*nTi[t]*nCO[t]) - (k20*nC[t]*nOH[t] + k21*nC[t]*nOH[t] + k22*nC[t]*nCO[t] + k26*nC[t]*nO2[t])), ..., nH[0] = 10^-7, nH2[0] = 10^-2, nC[0] = 10^-3, ... When solving this system, the program will give me the molecular concentration for each specie which I'd use as initial condition to the next iteration ; After solving the system, the program should go back to step 4, use the next value of t[r_,Φ_], calculate k1 to k143 again and use them to solve the system using the previously molecular concentration (nX) that was calculated before, until all values of t[r_,Φ_] are used.  Tsup[Φ_] := To + ΔT * Cos[Φ]; Rsup[Φ_] := Ro + ΔR * Cos[Φ + π]; To=2200 Ro=2.1 ΔT=440 ΔR=0.25  Having t[r_,Φ_] for each 'r' I can calculate 'k', wich goes from k1 to k143, to each t[r_,Φ_]: ki = α*(t[r,Φ]/300)^β * Exp[-ɣ/t[r,Φ]], where α,β and ɣ depend on the species. In the example above, these are the 'k' to solve the system of ODE: k1=(6.99*10^-14)*((t[r,Φ]/300)^2.8)*Exp[-1950/t[r,Φ]]; k2=2*10^-10 for all temperatures; k17=(3.14*10^-13)*((t[r,Φ]/300)^2.7)*Exp[-3150/t[r,Φ]]; k18=(2.05*10^-12)*((t[r,Φ]/300)^1.52)*Exp[-1736/t[r,Φ]]; k62=(1.77*10^-11)*Exp[178/t[r,Φ]]; k63=(1.85*10^-11)*((t[r,Φ]/300)^0.95)*Exp[-8571/t[r,Φ]]; k94 =(1.65*10^-12)*((t[r,Φ]/300)^1.14)*Exp[-50/t[r,Φ]]  Then I can finally calcute the molecular concentration profiles by solving a system of 50 non-linear, stiff, ordinary and coupled differential equations. Here an example of those equations.  s=NDSolve[{ nH'[t] == ((k17*nH2[t]*nO[t] + k18*nH2[t]*nOH[t] + k62*nO[t]*nOH[t]) - (k1*nH[t]*nOH[t] + k2*nH[t]*nH2O[t])), nH2'[t] == ((k1*nH[t]*nOH[t] + k2*nH[t]*nH2O[t]) - (k17*nH2[t]*nO[t] + k18*nH2[t]*nOH[t])), nO'[t] == ((k1*nH[t]*nOH[t] + k94*nOH[t]*nOH[t]) - (k62*nO[t]*nOH[t] + k63*nO[t]*nH2O[t])), nOH'[t] == ((k2*nH[t]*nH2O[t] + k17*nH2[t]*nO[t] + k63*nO[t]*nH2O[t]) - (k1*nH[t]*nOH[t] + k18*nH2[t]*nOH[t] + k62*nO[t]*nOH[t] + k94*nOH[t]*nOH[t])), nH2O'[t] == ((k18*nH2[t]*nOH[t] + k94*nOH[t]*nOH[t]) - (k2*nH[t]*nH2O[t] + k63*nO[t]*nH2O[t])), ..., nH[0] == 10^-7, nH2[0] == 5*10^-1, nO[0] == 10^-5, nOH[0] == 4.9*10^-7, nH2O[0] == 4*10^-6, }, ] When solving this system, the program will give me the molecular concentration for each specie which I'd use as initial condition to the next iteration ; After solving the system, the program should go back to step 4, use the next value of t[r_,Φ_], calculate k1 to k143 again and use them to solve the system using the previously molecular concentration (nX) that was calculated before, until all values of t[r_,Φ_] are used. 3 added 6 characters in body edited May 9 '15 at 2:50 André Oliveira 5355 bronze badges I believe I have a difficult problem to work on specially because I don't have any previously experience using Mathematica in this kind of situation. So I'll try to describe it as best as I can hoping that you can tell how should I proceed. Φ -> is the phase and I can freely set its initial value from 0 to 2π ; For each value of Φ, the program gives back a value to Tsup[Φ_] and Rsup[Φ_] :  Tsup[Φ_] := To + ΔT * Cos[Φ]; Rsup[Φ_] := Ro + ΔR * Cos[Φ + π];   Then I can calcute t[r_,Φ_] for each 'r' using both Tsup[Φ_] and Rsup[Φ_]. I can choose r's initial and maximum value, as well as the step: t[r_,Φ_] := Tsup[Φ]*[(1 - Sqrt[1 - (Rsup[Φ]/r)^2])]^4;]^1/4; Having t[r_,Φ_] for each 'r', I can calculate 'k', wich goes from k1 to k143, to each t[r_,Φ_]: ki = α*(t[r,Φ]/300)^βExp[-ɣt[r,Φ]] Then I can finally calcute the molecular concentration profiles by solving a system of 50 non-linear, stiff, ordinary and coupled differential equations. Here an example with three of those equations. nH'[t] == ((k17*nH2[t]*nO[t] + k18*nH2[t]*nOH[t] + k21*nC[t]*nOH[t] + k62*nO[t]*nOH[t] + k97*nOH[t]*nCO[t] + k141*nTi[t]*nOH[t]) - (k1*nH[t]*nOH[t] + k2*nH[t]*nH2O[t] + k5*nH[t]*nCO[t] + k8*nH[t]*nO2[t] + k10*nH[t]*nCO2[t])), nH2'[t] == ((k1*nH[t]*nOH[t] + k2*nH[t]*nH2O[t] + k143*nTi[t]*nH2O[t] + k144*nTiO[t]*nH2O[t]) - (k17*nH2[t]*nO[t] + k18*nH2[t]*nOH[t] + k19*nH2[t]*nO2[t])), nC'[t] == ((k5*nH[t]*nCO[t] + k140*nTi[t]*nCO[t]) - (k20*nC[t]*nOH[t] + k21*nC[t]*nOH[t] + k22*nC[t]*nCO[t] + k26*nC[t]*nO2[t])), ..., nH[0] = 10^-7, nH2[0] = 10^-2, nC[0] = 10^-3, ... When solving this system, the program will give me the molecular concentration for each specie which I'd use as initial condition to the next iteration ; After solving the system, the program should go back to step 4, use the next value of t[r_,Φ_], calculate k1 to k143 again and use them to solve the system using the previously molecular concentration (nX) that was calculated before, until all values of t[r_,Φ_] are used. I don't know if I made myself clear but I think it sums up the problem I have to solve. Can anyone help me by giving some tips and advices about how I should build the code? Of course that I don't want anyone to solve it for me. I'm learning how to work with Mathematica, so it would be very enlightening if someone help me to start to solve this problem. Thank you! I believe I have a difficult problem to work on specially because I don't have any previously experience using Mathematica in this kind of situation. So I'll try to describe it as best as I can hoping that you can tell how should I proceed. Φ -> is the phase and I can freely set its initial value from 0 to 2π ; For each value of Φ, the program gives back a value to Tsup[Φ_] and Rsup[Φ_] :  Tsup[Φ_] := To + ΔT * Cos[Φ]; Rsup[Φ_] := Ro + ΔR * Cos[Φ + π];  Then I can calcute t[r_,Φ_] for each 'r' using both Tsup[Φ_] and Rsup[Φ_]. I can choose r's initial and maximum value, as well as the step: t[r_,Φ_] := Tsup[Φ]*[(1 - Sqrt[1 - (Rsup[Φ]/r)^2])]^4; Having t[r_,Φ_] for each 'r', I can calculate 'k', wich goes from k1 to k143, to each t[r_,Φ_]: ki = α*(t[r,Φ]/300)^βExp[-ɣt[r,Φ]] Then I can finally calcute the molecular concentration profiles by solving a system of 50 non-linear, stiff, ordinary and coupled differential equations. Here an example with three of those equations. nH'[t] == ((k17*nH2[t]*nO[t] + k18*nH2[t]*nOH[t] + k21*nC[t]*nOH[t] + k62*nO[t]*nOH[t] + k97*nOH[t]*nCO[t] + k141*nTi[t]*nOH[t]) - (k1*nH[t]*nOH[t] + k2*nH[t]*nH2O[t] + k5*nH[t]*nCO[t] + k8*nH[t]*nO2[t] + k10*nH[t]*nCO2[t])), nH2'[t] == ((k1*nH[t]*nOH[t] + k2*nH[t]*nH2O[t] + k143*nTi[t]*nH2O[t] + k144*nTiO[t]*nH2O[t]) - (k17*nH2[t]*nO[t] + k18*nH2[t]*nOH[t] + k19*nH2[t]*nO2[t])), nC'[t] == ((k5*nH[t]*nCO[t] + k140*nTi[t]*nCO[t]) - (k20*nC[t]*nOH[t] + k21*nC[t]*nOH[t] + k22*nC[t]*nCO[t] + k26*nC[t]*nO2[t])), ..., nH[0] = 10^-7, nH2[0] = 10^-2, nC[0] = 10^-3, ... When solving this system, the program will give me the molecular concentration for each specie which I'd use as initial condition to the next iteration ; After solving the system, the program should go back to step 4, use the next value of t[r_,Φ_], calculate k1 to k143 again and use them to solve the system using the previously molecular concentration (nX) that was calculated before, until all values of t[r_,Φ_] are used. I don't know if I made myself clear but I think it sums up the problem I have to solve. Can anyone help me by giving some tips and advices about how I should build the code? Of course that I don't want anyone to solve it for me. I'm learning how to work with Mathematica, so it would be very enlightening if someone help me to start to solve this problem. Thank you! I believe I have a difficult problem to work on specially because I don't have any previously experience using Mathematica in this kind of situation. So I'll try to describe it as best as I can hoping that you can tell how should I proceed. Φ -> is the phase and I can freely set its initial value from 0 to 2π ; For each value of Φ, the program gives back a value to Tsup[Φ_] and Rsup[Φ_] :  Tsup[Φ_] := To + ΔT * Cos[Φ]; Rsup[Φ_] := Ro + ΔR * Cos[Φ + π];   Then I can calcute t[r_,Φ_] for each 'r' using both Tsup[Φ_] and Rsup[Φ_]. I can choose r's initial and maximum value, as well as the step: t[r_,Φ_] := Tsup[Φ]*[(1 - Sqrt[1 - (Rsup[Φ]/r)^2])]^1/4; Having t[r_,Φ_] for each 'r', I can calculate 'k', wich goes from k1 to k143, to each t[r_,Φ_]: ki = α*(t[r,Φ]/300)^βExp[-ɣt[r,Φ]] Then I can finally calcute the molecular concentration profiles by solving a system of 50 non-linear, stiff, ordinary and coupled differential equations. Here an example with three of those equations. nH'[t] == ((k17*nH2[t]*nO[t] + k18*nH2[t]*nOH[t] + k21*nC[t]*nOH[t] + k62*nO[t]*nOH[t] + k97*nOH[t]*nCO[t] + k141*nTi[t]*nOH[t]) - (k1*nH[t]*nOH[t] + k2*nH[t]*nH2O[t] + k5*nH[t]*nCO[t] + k8*nH[t]*nO2[t] + k10*nH[t]*nCO2[t])), nH2'[t] == ((k1*nH[t]*nOH[t] + k2*nH[t]*nH2O[t] + k143*nTi[t]*nH2O[t] + k144*nTiO[t]*nH2O[t]) - (k17*nH2[t]*nO[t] + k18*nH2[t]*nOH[t] + k19*nH2[t]*nO2[t])), nC'[t] == ((k5*nH[t]*nCO[t] + k140*nTi[t]*nCO[t]) - (k20*nC[t]*nOH[t] + k21*nC[t]*nOH[t] + k22*nC[t]*nCO[t] + k26*nC[t]*nO2[t])), ..., nH[0] = 10^-7, nH2[0] = 10^-2, nC[0] = 10^-3, ... When solving this system, the program will give me the molecular concentration for each specie which I'd use as initial condition to the next iteration ; After solving the system, the program should go back to step 4, use the next value of t[r_,Φ_], calculate k1 to k143 again and use them to solve the system using the previously molecular concentration (nX) that was calculated before, until all values of t[r_,Φ_] are used. I don't know if I made myself clear but I think it sums up the problem I have to solve. Can anyone help me by giving some tips and advices about how I should build the code? Of course that I don't want anyone to solve it for me. I'm learning how to work with Mathematica, so it would be very enlightening if someone help me to start to solve this problem. Thank you! 2 edited body edited May 8 '15 at 15:20 André Oliveira 5355 bronze badges I believe I have a difficult problem to work on specially because I don't have any previously experience using Mathematica in this kind of situation. So I'll try to describe it as best as I can hoping that you can tell how should I proceed. Φ -> is the phase and I can freely set its initial value from 0 to 2π ; For each value of πΦ, the program gives back a value to Tsup[Φ_] and Rsup[Φ_] :  Tsup[Φ_] := To + DeltaTCos[Φ_]; Rsup[Φ_] := Ro + DeltaRΔT * Cos[Φ]; Rsup[Φ_] := Ro + ΔR * Cos[Φ + π];  Then I can calcute T[r_t[r_,Φ_] for each 'r' using both Tsup[Φ_] and Rsup[Φ_]. I can choose r's initial and maximum value, as well as the step: t[r_,Φ_] := Tsup[Φ]*[(1 - Sqrt[1 - (Rsup[Φ]/r)^2])]^4; Having T[r_t[r_,Φ_] for each 'r', I can calculate 'k', wich goes from k1 to k143, to each T[r_t[r_,Φ_]: ki = α*(T[r_t[r,Φ_]Φ]/300)^βExp[-ɣT[r_t[r,Φ_]]Φ]] Then I can finally calcute the molecular concentration profiles by solving a system of 50 non-linear, stiff, ordinary and coupled differential equations. Here an example with three of those equations. nH'[t] == ((k17*nH2[t]*nO[t] + k18*nH2[t]*nOH[t] + k21*nC[t]*nOH[t] + k62*nO[t]*nOH[t] + k97*nOH[t]*nCO[t] + k141*nTi[t]*nOH[t]) - (k1*nH[t]*nOH[t] + k2*nH[t]*nH2O[t] + k5*nH[t]*nCO[t] + k8*nH[t]*nO2[t] + k10*nH[t]*nCO2[t])), nH2'[t] == ((k1*nH[t]*nOH[t] + k2*nH[t]*nH2O[t] + k143*nTi[t]*nH2O[t] + k144*nTiO[t]*nH2O[t]) - (k17*nH2[t]*nO[t] + k18*nH2[t]*nOH[t] + k19*nH2[t]*nO2[t])), nC'[t] == ((k5*nH[t]*nCO[t] + k140*nTi[t]*nCO[t]) - (k20*nC[t]*nOH[t] + k21*nC[t]*nOH[t] + k22*nC[t]*nCO[t] + k26*nC[t]*nO2[t])), ..., nH[0] = 10^-7, nH2[0] = 10^-2, nC[0] = 10^-3, ... When solving this system, the program will give me the molecular concentration for each specie which I'd use as initial condition to the next iteration ; After solving the system, the program should go back to step 4, use the next value of T[r_t[r_,Φ_], calculate k1 to k143 again and use them to solve the system using the previously molecular concentration (nX) that was calculated before, until all values of T[r_t[r_,Φ_] are used. I don't know if I made myself clear but I think it sums up the problem I have to solve. Can anyone help me by giving some tips and advices about how I should build the code? Of course that I don't want anyone to solve it for me. I'm learning how to work with Mathematica, so it would be very enlightening if someone help me to start to solve this problem. Thank you! I believe I have a difficult problem to work on specially because I don't have any previously experience using Mathematica in this kind of situation. So I'll try to describe it as best as I can hoping that you can tell how should I proceed. Φ -> is the phase and I can freely set its initial value from 0 to 2π ; For each value of π, the program gives back a value to Tsup[Φ_] and Rsup[Φ_] :  Tsup[Φ_] := To + DeltaTCos[Φ_]; Rsup[Φ_] := Ro + DeltaRCos[Φ + π];  Then I can calcute T[r_,Φ_] for each 'r' using both Tsup[Φ_] and Rsup[Φ_]. I can choose r's initial and maximum value, as well as the step: t[r_,Φ_] := Tsup[Φ]*[(1 - Sqrt[1 - (Rsup[Φ]/r)^2])]^4; Having T[r_,Φ_] for each 'r', I can calculate 'k', wich goes from k1 to k143, to each T[r_,Φ_]: ki = α*(T[r_,Φ_]/300)^βExp[-ɣT[r_,Φ_]] Then I can finally calcute the molecular concentration profiles by solving a system of 50 non-linear, stiff, ordinary and coupled differential equations. Here an example with three of those equations. nH'[t] == ((k17*nH2[t]*nO[t] + k18*nH2[t]*nOH[t] + k21*nC[t]*nOH[t] + k62*nO[t]*nOH[t] + k97*nOH[t]*nCO[t] + k141*nTi[t]*nOH[t]) - (k1*nH[t]*nOH[t] + k2*nH[t]*nH2O[t] + k5*nH[t]*nCO[t] + k8*nH[t]*nO2[t] + k10*nH[t]*nCO2[t])), nH2'[t] == ((k1*nH[t]*nOH[t] + k2*nH[t]*nH2O[t] + k143*nTi[t]*nH2O[t] + k144*nTiO[t]*nH2O[t]) - (k17*nH2[t]*nO[t] + k18*nH2[t]*nOH[t] + k19*nH2[t]*nO2[t])), nC'[t] == ((k5*nH[t]*nCO[t] + k140*nTi[t]*nCO[t]) - (k20*nC[t]*nOH[t] + k21*nC[t]*nOH[t] + k22*nC[t]*nCO[t] + k26*nC[t]*nO2[t])), ..., nH[0] = 10^-7, nH2[0] = 10^-2, nC[0] = 10^-3, ... When solving this system, the program will give me the molecular concentration for each specie which I'd use as initial condition to the next iteration ; After solving the system, the program should go back to step 4, use the next value of T[r_,Φ_], calculate k1 to k143 again and use them to solve the system using the previously molecular concentration (nX) that was calculated before, until all values of T[r_,Φ_] are used. I don't know if I made myself clear but I think it sums up the problem I have to solve. Can anyone help me by giving some tips and advices about how I should build the code? Of course that I don't want anyone to solve it for me. I'm learning how to work with Mathematica, so it would be very enlightening if someone help me to start to solve this problem. Thank you! I believe I have a difficult problem to work on specially because I don't have any previously experience using Mathematica in this kind of situation. So I'll try to describe it as best as I can hoping that you can tell how should I proceed. Φ -> is the phase and I can freely set its initial value from 0 to 2π ; For each value of Φ, the program gives back a value to Tsup[Φ_] and Rsup[Φ_] :  Tsup[Φ_] := To + ΔT * Cos[Φ]; Rsup[Φ_] := Ro + ΔR * Cos[Φ + π];  Then I can calcute t[r_,Φ_] for each 'r' using both Tsup[Φ_] and Rsup[Φ_]. I can choose r's initial and maximum value, as well as the step: t[r_,Φ_] := Tsup[Φ]*[(1 - Sqrt[1 - (Rsup[Φ]/r)^2])]^4; Having t[r_,Φ_] for each 'r', I can calculate 'k', wich goes from k1 to k143, to each t[r_,Φ_]: ki = α*(t[r,Φ]/300)^βExp[-ɣt[r,Φ]] Then I can finally calcute the molecular concentration profiles by solving a system of 50 non-linear, stiff, ordinary and coupled differential equations. Here an example with three of those equations. nH'[t] == ((k17*nH2[t]*nO[t] + k18*nH2[t]*nOH[t] + k21*nC[t]*nOH[t] + k62*nO[t]*nOH[t] + k97*nOH[t]*nCO[t] + k141*nTi[t]*nOH[t]) - (k1*nH[t]*nOH[t] + k2*nH[t]*nH2O[t] + k5*nH[t]*nCO[t] + k8*nH[t]*nO2[t] + k10*nH[t]*nCO2[t])), nH2'[t] == ((k1*nH[t]*nOH[t] + k2*nH[t]*nH2O[t] + k143*nTi[t]*nH2O[t] + k144*nTiO[t]*nH2O[t]) - (k17*nH2[t]*nO[t] + k18*nH2[t]*nOH[t] + k19*nH2[t]*nO2[t])), nC'[t] == ((k5*nH[t]*nCO[t] + k140*nTi[t]*nCO[t]) - (k20*nC[t]*nOH[t] + k21*nC[t]*nOH[t] + k22*nC[t]*nCO[t] + k26*nC[t]*nO2[t])), ..., nH[0] = 10^-7, nH2[0] = 10^-2, nC[0] = 10^-3, ... When solving this system, the program will give me the molecular concentration for each specie which I'd use as initial condition to the next iteration ; After solving the system, the program should go back to step 4, use the next value of t[r_,Φ_], calculate k1 to k143 again and use them to solve the system using the previously molecular concentration (nX) that was calculated before, until all values of t[r_,Φ_] are used. I don't know if I made myself clear but I think it sums up the problem I have to solve. Can anyone help me by giving some tips and advices about how I should build the code? Of course that I don't want anyone to solve it for me. I'm learning how to work with Mathematica, so it would be very enlightening if someone help me to start to solve this problem. Thank you! 1 asked May 8 '15 at 14:14 André Oliveira 5355 bronze badges