Plotting a table with many y values

Question

I'm trying to plot a bifurcation map for a Lorenz System. I have a function (called estables in the code that follows) which, for a certain value r, returns the list of stationary points for the system with that parameter r.

Then I do this:

r0 = 0;
rf = 10;
dr = 0.5;
ys = Table[estables[ra], {ra, r0, rf, dr}]

ListPlot[ Flatten[ys, {2}], PlotRange -> All , DataRange -> {r0, rf}, 
PlotStyle -> {Blue, PointSize[Small]}, Joined -> False ]

And I get a weird thing where the points are plotted out of order. That is, the points in the first line are not the ones plotted above the 0, as it should be.

What would be the best way to achieve this?

Edit: This are the differential equations:

ecf = σ (p[t] - f[t]);
ecp = -p[t] + d[t] f[t];
ecd = b (r - d[t] - f[t] p[t]);

and this is the function

t0 = 100;
tf = 105;
dt = 0.05;

estables[ra_] :=

Module[{par, solnum, fs}, 
    par = {σ -> 3., b -> 1, r -> ra}; 
    solnum = 
      NDSolve[{Derivative[1][f][t] == ecf, 
               Derivative[1][p][t] == ecp, 
               Derivative[1][d][t] == ecd, 
               f[0] == 0.00001, p[0] == 0., 
               d[0] == 1} /. par, 
               {f, p, d}, {t, t0, tf}, 
               MaxSteps -> 100000
    ];
    fs = Flatten[Table[Evaluate[f[t] /. solnum], {t, t0, tf, dt}]]; 
    Union[FindPeaks[fs], -FindPeaks[-fs]]][[;; , 2]]

Answer 1

Thanks to J.M. I found a more elegant way to obtain the statinoary values that returned far simpler data, which was easy to plot. This is how my code looks now:

ecf = s (p[t] - f[t]);
ecp = -p[t] + d[t] f[t];
ecd = b (r - d[t] - f[t] p[t]);

par = {s -> 3., b -> 1, r -> 100};
t0 = 800;
tf = 850;



estables[ra_] :=

 Module[{par, solnum,  points}, 
  par = {s -> 3., b -> 1, r -> ra}; 
  {solnum, points} = Reap@NDSolve[{
           Derivative[1][f][t] == ecf, 
           Derivative[1][p][t] == ecp, 
           Derivative[1][d][t] == ecd, 
           f[0] == 0.00001, p[0] == 0., 
           d[0] == 1, 
           WhenEvent[{f'[t] == 0 && t > t0}, Sow[f[t]]]} 
           /. par, {f, p, d}, {t, t0, tf}, MaxSteps -> 100000];
  Flatten[points]
];


r0 = 0;
rf = 105;
dr = 0.2;

ys = Table[{ra, e}, {ra, r0, rf, dr}, {e, estables[ra]}];
ListPlot[ Flatten[ys, 1], PlotRange -> All , DataRange -> {r0, rf}, 
 PlotStyle -> {Blue, PointSize[0.002]}, Joined -> False, 
 AxesLabel -> {r, f} ]