# Bifurcation diagram for Lorenz attractor

As I was following a similar question I came across a beautiful answer here. However, in that code, he used InternalBag[] which heavily depends on ram that I can't afford and I need to make a data file too. So Please Help me edit this code so that I can export the data to my hard disk.

Along with that, this code only stores the extremum values of Z as in WhenEvent there is only one condition i.e. z'[t]==0. I want to store minimum and maximum values in two different columns. How do I do that?

res = InternalBag[];(*a place to store results*)
tmax = 100;(*how long to run for each r value*)
{x0, y0, z0} = {10, 10, 10};(*initial ICs*)

Do[sol = With[{a = 10.0, c = 8/3},
NDSolve[{x'[t] == a (y[t] - x[t]),
y'[t] == r x[t] - y[t] - x[t] z[t], z'[t] == x[t] y[t] - c*z[t],
x[0] == x0, y[0] == y0, z[0] == z0,(*save extrema of z[t]*)
WhenEvent[z'[t] == 0, InternalStuffBag[res, {r, z[t]}]]}, {x,
y, z}, {t, 0, tmax}, Method -> "StiffnessSwitching",MaxSteps -> \[Infinity]]][[1]];
(*save end value for next ICs*){x0, y0, z0} = {x[tmax], y[tmax], z[tmax]} /. sol;, {r, 100, 10, -0.1}];

ListPlot[InternalBagPart[res, All], PlotStyle -> {Red, Opacity[0.1], PointSize[0.001]}, AxesLabel -> {"r", "z"}]


I was suggested to use Reap or Sow, but I am unable to edit the code properly

• Glad you got an answer to your question. BTW, where did you see that InternalBag is obsolete? Aug 23, 2022 at 12:26
• @ChrisK Sorry, I was wrong there. It's not obsolete but it heavily depends on ram. Aug 24, 2022 at 3:40

Here is an example how you would use Sow and Reap. Note, I decreased the r interval, as I did not want to wait a long time.

tmax = 100;(*how long to run for each r value*)
{x0, y0, z0} = {10, 10, 10};(*initial ICs*)

res = Reap[
Do[sol =
With[{a = 10.0, c = 8/3},
NDSolve[{x'[t] == a (y[t] - x[t]),
y'[t] == r x[t] - y[t] - x[t] z[t],
z'[t] == x[t] y[t] - c*z[t], x[0] == x0, y[0] == y0,
z[0] == z0,(*save extrema of z[t]*)
WhenEvent[z'[t] == 0, Sow[{r, z[t]}]]}, {x, y, z}, {t, 0,
tmax}, Method -> "StiffnessSwitching",
MaxSteps -> \[Infinity]]][[1]];
(*save end value for next ICs*){x0, y0,
z0} = {x[tmax], y[tmax], z[tmax]} /. sol;, {r, 11,
10, -0.1}]][[2]];

ListPlot[res, PlotStyle -> {Red, Opacity[0.1], PointSize[0.001]},
AxesLabel -> {"r", "z"}]

• Thank you. How do I separate minimum and maximum values in WhenEvent? Aug 23, 2022 at 9:15
• In the WhenEvent you only have access to 3 data: t, z[t] and z'[t]. That is not enough to determine if we have a min or max. But "sol" contains the whole function and you may use it to determine if it is min or max. Aug 23, 2022 at 10:00