This notebook finds an approximate initial condition that is in the stable manifold of the origin.  Of course the whole z-axis in in the stable manifold of the origin. 
But we want to find a point in the Lorenz attractor that is in the stable manifold of the origin.
We have a target of 1/2, which has the binary sequence {1, 0}.  (It goes to the right once and then goes straight into the origin.
The two observable solutions have sign sequences
{1,-1, 1, 1, 1 (repeated 1s.) } and
{1,1,-1,-1, (repeated -1s.) }

The first cell integrates the Lorenz equation for 0 < t < tf, with the initial conditions on the line segment jointing the two nontrivial fixed points.
The parameter s, between -1 and 1, determines the position on the line segment.

sig = 10;
r = 28;
b = 8/3;
xf = Sqrt[(r-1.) b];
zf = r-1;

(* The fixed points are (0,0,0), (xf, xf, zf), and (-xf, -xf, zf) *)
fig[s_, tf_] := Module[{},
x0 = s* xf;  
y0 = s* xf;
z0 = r-1;
{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[z[t] > 27, Sow[{x[t],y[t], t}] ]},
{x,y, z},{t,0,tf}];
soln = Part[soln, 1];
If[poincare ≠ {}, poincare = Part[poincare, 1] ];
points = Table[{poincare[[i, 1]], poincare[[i, 2]], r-1}, {i, Length[poincare]}];
yList = Table[poincare[[i, 2]], {i, Length[poincare]}];
signList = Sign[yList];
sum = 0.;
Do[sum = sum + Sign[poincare[[i, 1]] ] /2^(i), {i, 1, Length[poincare]}];
z27plane = Graphics3D[{Opacity[0.7],Polygon[{{-30, -30, r-1}, {-30, 30, r-1}, {30, 30, r-1}, {30, -30, r-1}}], PointSize[Medium], Black, Point[{xf, xf, r-1}], Point[{-xf, -xf, r-1}],Point[{x0, y0, z0}],Red, Point[points]}];
pTrajectory =ParametricPlot3D[{x[t], y[t], z[t]}/.soln, {t, 0, tf}, PlotLabel -> "IC = (s (Subscript[x, f]), s (Subscript[x, f]), 27), with s = "<>ToString[s]<>"\nsignList = "<>ToString[signList]<>"\nf(s) = "<>ToString[sum]];
Show[pTrajectory, z27plane, PlotRange -> All]
]

The next two cells show two very close initial conditions, one which turns to our left, and the other which turns to our right.

fig[0.17089041,10]

fig[0.17089046,10]

poincare (* List of {x,y,t} when z = 27, with z increasing. *)

The following function just returns f(s) (for a finite time run 0 < t < tf) and does not do the graphics.

f[s_, tf_] := Module[{},
If[s == -1, Return[-1] ];
If[s == 0, Return[0] ];
If[s == 1, Return[1]];
x0 = y0 =s xf;
z0 = r-1;
{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[z[t] > 27, Sow[{x[t],y[t], t}] ]},{x,y, z},{t,tf}];
soln = Part[soln, 1];
poincare = Part[poincare, 1];
signList = Sign[Transpose[poincare ][[1]] ];
sum = 0.;
Do[sum = sum + signList[[i]] /2^(i), {i, 1, Length[poincare]}];
sum
]

f[.1, 10]
signList

poincare

The next cell shows the graph of the function f.  Note how quickly the function saturates to 1 or -1.

pts = {};
Do[AppendTo[pts, {s, f[s, 10]}], {s, -1, 1, .02}]
ListPlot[pts]

I can use this to (hopefully) find orbits with any sequence I want.  I can just use root-finding.  Choose a sequence like LRLLR LRLLR LRLLR ... and get the number for that, and Use FindRoot with an interval initial guess.

s1 = 0.17;
s2 = .172;
tf = 18;
maxIts = 20;
tol = 0.0000001;
target = 1/2;
f1 = f[s1, tf] - target
Print[signList]
f2 = f[s2, tf]-target 
Print[signList];
its = 0;
While[Abs[f2 -f1] > tol && its < maxIts,
snew = s1 - f1*(s2-s1)/(f2-f1);
fnew = f[snew, tf] - target;
Print["fnew = ",fnew];
Print["snew = ",snew];
Print[signList];
s2 = s1;
f2 = f1;
s1 = snew;
f1 = fnew;
its++
]

Here are the last two s values, at “full resolution”.  The true value of s  = s*is somewhere between them.  With that value of s, the solution goes to the origin

0.17089041133604407`
0.170890457637069`

x0 = snew*xf;
y0 = snew*xf;
Print["The initial condion value of s that approaches origin after going right is "]
snew
z0 = r-1;
soln = 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},{x,y, z},{t,0, tf}];

 ParametricPlot3D[{x[t], y[t], z[t]}/.soln, {t, 0, tf}, PlotRange -> All]

x0 = snew*xf;
y0 = snew*xf;
Print["The initial condion value of s that approaches origin after going right is ", snew]
snew
z0 = r-1;
tPlot = 5;
soln = 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},{x,y, z},{t,0,tPlot}];

 ParametricPlot3D[{x[t], y[t], z[t]}/.soln, {t, 0, tPlot}, PlotRange -> All]