# WhenEvent, Storing Sow data, Loops

My problem is essentially to do with manipulating and storing Sow data in a nested loop structure.

I have a set of differential equations describing time-evolution of a bunch of variables. Eventually I want the initial conditions i0, j0, k0 to be parameters. I have a range for each of these that I want to scan over. Eventually, I want to find out the t value at which point i[t]==1.

So far I can set up an NDSolveValue taking parameteres, and Reap for a specific set of parameters.

test[i0_, j0_, k0_] := NDSolveValue[{coupled nonlinear first order ODEs,
i[0] == i0, j[0] == j0, k[0] == k0, WhenEvent[i[t] == 1, Sow@{i0, j0, k0};
Sow[t]]}, {i, j, k}, {t, 0, 15}];

test[0.4, 0.001, 0.03]; // Reap


The output makes sense, giving the t-value as 11.2293 for this case which seems correct:

{Null, {{{0.4, 0.001, 0.03}, 11.2293}}}


I have been struggling with scanning across a specified i0, j0, k0 range. Here is something that kind of works (StopIntegration is to avoid a pole that happens at some late t):

Table[Table[Table[
NDSolveValue[{the ODEs, i[0] == i0, j[0] == j0, k[0] == k0,
WhenEvent[i[t] == 1, Print[t,i0,j0,k0];
Sow[t, {i0,j0,k0}]; "StopIntegration"]}, {i, j, k}, {t, 0, 25}],
{i0,0.3, 0.4, 0.1}], {j0, 0.0001, 0.005, 0.001}], {k0, 0.1, 0.5, 0.1}]


And this gives me:

22.1212....0.3....0.0001....0.1

11.2014....0.4....0.0001....0.1

6.79366....0.5....0.0001....0.1

22.1209....0.3....0.0021....0.1

for each combination of values, which also seem correct. However, I am struggling with saving these results in a proper data (.dat) file. I want to ultimately make a table with columns i0, j0, j0, t which can be exported as a .dat file. Ideally, I would also like to append to this table the columns f(i0,j0,k0), g(i0,j0,k0), which are certain functions of the values. How do I create a professional data table out of this data (I've been using 'Print' because I don't know how to use Sow/Reap data properly)?

Edit: If I add Reap@ to the start of the expression it gives me a list of t-values for i0=i1 then for i0=i2 then i0=i3, that is, for each constant i from the i range, with the others varying; followed by a list of t-values for j0=j1,....

At this stage my query is mainly to do with formatting the data, i.e. generating a neat exportable table. How do I convert the Reap output into a data table?

Any help would be appreciated!

• Can you be more concrete and specific? ijk={{.4,.001,.03},{.6,.002,.06},{.8,.004,.1}} Now run my code once for each of those three triples? Or run my code for every possible subset of those 9 values? Can you show precisely how you want the result to be a matrix of values from the 3 runs? Maybe skip getting confused with writing the code and just clearly tell someone exactly what calculations you want them to do and exactly how you want them to save the results.
– Bill
Jun 14, 2018 at 18:07
• @Bill - I have added some more detail. Using nested Tables, I am able to loop across the i0, j0, and k0 ranges that I want (see edited code). The Print function seems to verify that it is working. Using Reap@ at the start of the expression I get lists of the form { {tlist(i1)}, {tlist(j1), {tlist(k1)}, {tlist(i2)}, {tlist (j2)}...} where by tlist(i1) I mean the t-value for all the permutations with i0=i1 (first element in i-list). This has all the data I need, except I would ultimately like a table with columns t, i0, j0, k0. Jun 14, 2018 at 18:32
• @Bill - As for your question, I would ideally like to know how to make the code work with every 'subset'/permutation as in the current range method, although, it would also be useful to eventually have a way to do it for a list of triplets, as that would be useful for the physical system I am considering. As long as ultimately I am still able to generate that kind of data table that I mentioned. Thank you very, very much for your help! Jun 14, 2018 at 18:34
• Change your Reap[Table[Table[Table[stuff;Sow[{t,{i,j,k}}],{i,..}],{j,..}],{k,..}]] to nBy4Matrix=Reap[Do[stuff;Sow[{t,i,j,k}],{i,...},{j,...},{k,...}]][[2,1]] and see if you get a nice n by 4 matrix for all your i,j,k values. Then you can do MatrixForm[nBy4Matrix] and see if that gives you your formatting.
– Bill
Jun 14, 2018 at 23:14
• @Bill - Thanks, that actually solved my problems! Would you like to compile it into an answer so that I may accept and upvote it? Also, if I had a list of {{i1, j1, k1}, {i2, j2, k2}....} values, how would I pass that? Any words to illuminate the syntax choice and Do over Table would also really help :) Jun 15, 2018 at 0:57

There are always at least a dozen different ways of doing anything in Mathematica, and that doesn't even count turning function names into punctuation. Here are three ways.

Assume that stuff does your calculations, ends up assigning a value to t and that multiple steps are separated by semicolons.

Table with multiple iterators will create nested lists, but Flatten with an appropriate 2nd argument will reduce that to a simple matrix.

nBy4Matrix = Flatten[Table[stuff; {t, i, j, k}, {i,1,3}, {j,1,3}, {k,1,3}], 2]


Sow and Reap will let you accumulate results from any sequence of calculations.

nBy4Matrix = Reap[Do[stuff; Sow[{t, i, j, k}], {i,1,3}, {j,1,3}, {k,1,3}]][[2,1]]


Map will let you perform a calculation on each item from a list, each of which happens to be a triple in your case, and return a list of the results to you.

nBy4Matrix = Map[({i,j,k}=#; stuff; {t, i, j, k})&, {{1,1,1}, {2,3,1}, {3,1,3}, {2,1,1}}]


Pick one way that you will be able to remember in the future and be able to use without making any mistakes. That is the real goal of all this.

• This may be a silly question, but -- what if I wanted some function of i/j/k? For example, Sow[{t,4*i,Sqrt[j],i*k}] performs the operation and gives me what I want -- but my question is, what if I want to scan across equally spaced steps in 4i, Sqrt[j], etc? Simply modifying them as so while specifying the i,j,k range doesn't work.,, Jun 23, 2018 at 15:22