Here's an example exhibited by Nusc, which I have tweaked slightly to make it look more like your mathematica example.

### Reference: http://www.mapleprimes.com/questions/36580-Bifurcation-Diagram

### xexpr is the logistic function to be iterated (we always start off at x=1/2, which will eventually attract).

### [ra,rb] is the range of the parameter.

### acc is the number of points sampled in [ra,rb]

Bifurcation := proc(initialpoint,xexpr,ra,rb,acc)

local p1,hr,A,L1,i,j,phi:

global r,L2:

hr := unapply(xexpr,x);

A := Vector(600):

L1 := Vector(acc*500):

for j from 1 to acc+1 do

r := (ra + (j-1)*(rb-ra)/acc):

A[1] := hr(initialpoint):

for i from 2 to 500 do

A[i] := evalf(hr(A[i-1])):

end do:

for i from 1 to 400 do

L1[i+400*(j-1)] := [r,A[i+100]]:

end do:

end do:

L2 := {seq(L1[i], i = 1..acc*400)}:

p1 := plots:-pointplot(L2, 'symbol' = solidcircle, 'symbolsize' = 8, 'color' = blue):

unassign('r'):

return(p1):

end proc:

### Example

P1 := Bifurcation(1/2,r*x*(1-x),2.5,4,250):

plots:-display(P1, 'axes' = box, 'labels' = [r, x] );

And the second graph is from the wikipedia page, quite pretty:

http://en.wikipedia.org/wiki/File:LogisticMap_BifurcationDiagram.png