# Using Mathematica to solve a recursive system of differential equations [closed]

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.

1. Φ -> is the phase and I can freely set its initial value from 0 to 2π ;

2. For each value of Φ, the program gives back a value to Tsup[Φ_] and Rsup[Φ_] :

 Tsup[Φ_] := To + ΔT * Cos[Φ];
Rsup[Φ_] := Ro + ΔR * Cos[Φ + π];
To=2200
Ro=2.1
ΔT=440
ΔR=0.25


1. 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;

1. 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,Φ]] 

1. 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,
},
]

1. 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 ;

2. 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!

## closed as off-topic by user9660, m_goldberg, MarcoB, Yves Klett, ÖskåDec 24 '15 at 10:44

• The question does not concern the technical computing software Mathematica by Wolfram Research. Please see the help center to find out about the topics that can be asked here.
If this question can be reworded to fit the rules in the help center, please edit the question.

• Welcome to Mathematica.SE! I suggest that: 1) You take the introductory Tour now! 2) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! 3) As you receive help, try to give it too, by answering questions in your area of expertise. – bbgodfrey May 8 '15 at 14:27
• In step 2, you say π when you may mean Φ. Also, you use t[r_,Φ_]  in some places and T[r_,Φ_]  in others. Please make the necessary corrections. Finally, what are DeltaTCos and DeltaRCos. – bbgodfrey May 8 '15 at 14:37
• Thank you for the suggestions. I've made the necessary corrections. – André Oliveira May 8 '15 at 15:25
• You may obtain some responses, if you provide code that readers can attempt to run. That includes, provide (say) three ODEs that contain only three unknowns, and define all constants. Since stiffness probably is you primary issue, select equations that you know to be stiff. – bbgodfrey May 9 '15 at 14:47
• I'm voting to close this question as off-topic because it is too localized; i.e, it applies only to the local situation and needs of its poster and answers will not benefit others. – m_goldberg Dec 23 '15 at 6:29