The Lorenz equations are
dx/dt = sigma (Y - X)	
dy/dt =	-XZ + rX - Y	
dz/dt = XY - bZ
We consider the standard parameters, r = 28, sigma = 10, b = 8/3.
The fixed points are (0, 0, 0), (xf, xf, 27) and (-xf, -xf, 27), 

Two functions are defined when you run the next cell.
The inputs for each function are s (best chosen between 0 and 1) and tf (the final time of the numerical integration).

The initial conditions below are chosen on the line segment between the two fixed point on the z = 27 plane:
(s xf, s xf, 27)

The fig[s, tf] command shows the Poincare map, with a dot every time the trajectory crosses the z = 27 plane going upward.

The list[s, tf] command gives a list of 1’s (for looping right) and -1’s (for looping left).

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) *)

fig[tf_] := Module[{},
x0 = 0.046;  (* the intial condition near origin *)
y0 = 0.1;
z0 = 0;
{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](* Even though the sulution is unique, the NDSolve function gives a list of solutions. Take the first and only one*);
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];
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 near (0,0,0)"];
Show[pTrajectory, z27plane, PlotRange -> All]
]

list[tf_] := Module[{},
x0 = 0.046;  
y0 = 0.1;
z0 = 0;
{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]
]

Now that these functions are defined, you can run the functions. 
The parameters are
* s between -1 and 1 to give the initial condition
 * tf, or t-final.  The numerical integration is done for 0 < t < tf.  You can make tf quite large, although the figure looks best with moderate tf.
 
The figure shows the plane z = 27, which contains the two fixed points (± (Subscript[x, f]), ± (Subscript[x, f]), 27) which are shown as black dots. Another black dot is the initial condition which is on the line segment between the two fixed points for s between -1 and 1. The red dots show where the trajectory pierces the plane z = 27 while going upward.  The  fig[s, tf] function shows the upward piercings for 0 < t < tf.  The number of piercings is approximately tf.  The list[s, tf] function keeps track of the piercing happens at positive or negative x.

When s is close to -1 or 1, the initial condition is very close to a fixed point, and it spins around many times before crossing to the other side.  For example, with s = 0.9 it crosses after 20 upward hits on the right.
 After making that first crossing to the other side, it is essentially on the Lorenz attractor.

tf =18.5;
fig[tf]
list[tf]
list[100]

The sign sequence of the unstable manifold has 25 “-1”s in a row.

Length[{1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,1,1}]

Note :  If you are doing statistics on the list of L' s and R' s (-1 s and 1 s) with a random s between 0 and 1 (s = Random[]), then disregard the initial string of repeated 1 s .

Df = {{-sig, sig, 0},{r, -1, 0}, {0, 0, -b}};
Print["Df = ", MatrixForm[Df]]
Eigenvalues[Df]//N
{evals, evecs} =Eigensystem[Df]
{evals, evecs} =Eigensystem[Df]//N