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