This notebook makes a dense map of z_n+1 vs z_n.
Here z_n is the nth local maximum value of z(t) on a solution to the Lorenz equations.
We have found two values of s that just miss going to the origin.
These are two approximations of the solution with f(s) =1/2, Sequence {R,0} 
The first, sRm has f(sRm) = 1/2 - epsilon
and the second, sRp has f(sRp) = 1/2 + epsilon

sR0m = 0.170890443;
sR0p = 0.170890444;

sig = 10;
r = 28;
b = 8/3;
(* You can change these constants above this line. *)

xf = Sqrt[(r-1.) b];
(* The fixed points are (0,0,0), (xf, xf, r-1), and (-xf, -xf, r-1) *)

sToZPairs[s_,tf_] := Module[{zPairList},
x0 = s*xf;  
y0 =s*xf;
z0 = 27;
{soln,poincare}=Reap@NDSolve[{
x'[t]==sig(y[t]-x[t]),
y'[t]==r x[t]-y[t]- x[t] z[t], 
z'[t] == x[t] y[t] - b z[t],
x[0]== x0,y[0]==y0,z[0] == z0, 
WhenEvent[x[t]*y[t]-b*z[t] < 0, Sow[z[t]] ]},
{x,y, z},{t,0,tf}];
If[poincare ≠ {}, poincare = Part[poincare, 1] ];
zList = Table[poincare[[i]], {i, Length[poincare]}];
zPairList = Table[{zList[[i]], zList[[i+1]]}, {i, 1, Length [zList]-1}]
]

zMaxICPairList[IC_,tf_] := Module[{},
{x0,y0,z0} = IC;
{soln,poincare}=Reap@NDSolve[{
x'[t]==sig(y[t]-x[t]),
y'[t]==r x[t]-y[t]- x[t] z[t], 
z'[t] == x[t] y[t] - b z[t],
x[0]== x0,y[0]==y0,z[0] == z0, 
WhenEvent[x[t]*y[t]-b*z[t] < 0, Sow[z[t]] ]},
{x,y, z},{t,0,tf}];
If[poincare ≠ {}, poincare = Part[poincare, 1] ];
zList = Table[poincare[[i]], {i, Length[poincare]}];
zPairList = Table[{zList[[i]], zList[[i+1]]}, {i, 1, Length [zList]-1}]
]

zMaxICPairList[{.0046,.01,0}, 3]

zU = 48.36162248297621`; (* Copied from output of previous command *)

tf = 20;
zListRp =sToZPairs[sR0p, tf]
Print["(Subscript[z, u]) = ", zU, " the largest zMax, the first zMax on the unstable manifold."];
zS = zListRp[[1,1]];
Print["(Subscript[z, s]) = ", zS, " headed to the origin, the frist zMax on the stable manifold."];
zL  = zListRp[[2,2]];
Print["(Subscript[z, l]) = ", zL, " the zMax of the first loop on the unstable manifold."];
zU-zS
zS - zL
p1 =ParametricPlot3D[{x[t], y[t], z[t]}/.soln, {t, 0, tf}, PlotRange -> All]

mapPairs = zListRp;
tf = 5;
Do[zList = sToZPairs[s, tf];
mapPairs = Union[mapPairs, zList],
{s,.1,.8,.01 }
];
Do[zList = sToZPairs[s, tf];
mapPairs = Union[mapPairs, zList],
{s,sR0p-.1,sR0p+.1,.01 }];
Do[zList = sToZPairs[s, tf];
mapPairs = Union[mapPairs, zList],
{s,sR0p-.009,sR0p+.011,.001 }];tf = 10;Do[zList = sToZPairs[s, tf];
mapPairs = Union[mapPairs, zList],
{s,sR0p-.0009,sR0p+.0001,.0001 }];
Do[zList = sToZPairs[s, tf];
mapPairs = Union[mapPairs, zList],
{s,sR0p-.00009,sR0p+.00001,.00001 }];
Do[zList = sToZPairs[s, tf];
mapPairs = Union[mapPairs, zList],
{s,sR0p-.000009,sR0p+.000001,.000001 }];
p1 = Plot[x, {x, 27, 49}, AspectRatio -> Automatic, PlotStyle -> Thin, PlotLabel -> "(Subscript[z, n + 1]) vs (Subscript[z, n])",Prolog->{LightYellow,Rectangle[{zL, zL},{zU, zU}]},Ticks ->{{{zL, "(Subscript[z, l])"}, 30,35,{zS ,"(Subscript[z, s])"},40, 45,{zU, "(Subscript[z, u])"}, {zL, "(Subscript[z, l])"}},{{zL, "(Subscript[z, l])"}, 30,35,{zS ,"(Subscript[z, s])"},40, 45,{zU, "(Subscript[z, u])"}}, {27, {zL, "(Subscript[z, l])"}}}];
pZmap=ListPlot[mapPairs, PlotRange -> All, AspectRatio -> Automatic, Ticks ->{{{zL, "(Subscript[z, l])"}, 30,35,{zS ,"(Subscript[z, s])"},40, 45,{zU, "(Subscript[z, u])"}, {zL, "(Subscript[z, l])"}},{{zL, "(Subscript[z, l])"}, 30,35,{zS ,"zs"},40, 45,{zU, "zu"}}, {27, {zL, "(Subscript[z, l])"}}}];
LorenzMap =Show[p1, pZmap]

mapPairs = zListRp;
tf = 5;
Do[zList = sToZPairs[s, tf];
mapPairs = Union[mapPairs, zList],
{s,.1,.8,.01 }
];
Do[s= sR0p + .3(i/100)^5;
zList = sToZPairs[s, tf];
mapPairs = Union[mapPairs, zList],
{i,-100,100,1 }];
p1 = Plot[x, {x, 27, 49}, AspectRatio -> Automatic, PlotStyle -> Thin, PlotLabel -> "(Subscript[z, n + 1]) vs (Subscript[z, n])",Prolog->{LightYellow,Rectangle[{zL, zL},{zU, zU}]},Ticks ->{{{zL, "(Subscript[z, l])"}, 30,35,{zS ,"(Subscript[z, s])"},40, 45,{zU, "(Subscript[z, u])"}, {zL, "(Subscript[z, l])"}},{{zL, "(Subscript[z, l])"}, 30,35,{zS ,"(Subscript[z, s])"},40, 45,{zU, "(Subscript[z, u])"}}, {27, {zL, "(Subscript[z, l])"}}}];
pZmap=ListPlot[mapPairs, PlotRange -> All, PlotStyle ->PointSize[0.004],AspectRatio -> Automatic, Ticks ->{{{zL, "(Subscript[z, l])"}, 30,35,{zS ,"(Subscript[z, s])"},40, 45,{zU, "(Subscript[z, u])"}, {zL, "(Subscript[z, l])"}},{{zL, "(Subscript[z, l])"}, 30,35,{zS ,"zs"},40, 45,{zU, "zu"}}, {27, {zL, "(Subscript[z, l])"}}}];
LorenzMap =Show[p1, pZmap]

Export["/Users/james/Library/CloudStorage/Dropbox/classes/MAT 667 2025.1/LorenzMap.pdf", LorenzMap]