Render bifurcation diagram (octave script)

bifurcation.m

## Period doubling bifurcation
## https://en.wikipedia.org/wiki/Logistic_map

1;

function ret = popfn(prev, r)
  ret = r * (prev) * (1-prev);
endfunction

# Iterate till n-iterations
function ret = popfniter(numiter, r, initial)
  nv = ones(1, numiter+1) .* initial;
  for iter = 2:numiter+1
    nv(iter) = popfn(nv(iter-1), r);
    if nv(iter) > 10000,
      nv(iter) = 10000;
    elseif nv(iter) < -10000,
      nv(iter) = -10000;
    endif;
  endfor;
  ret = arrayfun(@(x) round(x*10000)/10000, nv);
endfunction

## List of converging values
function ret = popfconv(numiter, r, initial)
  ll = popfniter(numiter, r, initial);
  [u, ~, j] = unique(ll);
  aa = [accumarray(j', 1), u'];
  aa = sortrows(aa, 1);
  numelems = round((numiter+1)/aa(end,1))-1;
  ret = aa(end-numelems:end);
endfunction;

## list of points.
function ret = poplot(numiter, rrange, initial)
  points = [];
  for r = rrange
    m = popfconv(numiter, r, initial);
    points = [points; [ones(1, length(m)) * r; m]'];
  endfor
  ret = points;
endfunction;

m = poplot(2000, 1:0.02:4, 0.5);
scatter(m(:,1)', m(:,2)', '.')