r_min = 2.8; r_max = 4; # the range of parameter values we study
n = 500; # the number of parameter values we consider in this range

t_max = 500; # how many iterations to simulate per parameter value
p_max = 100; # the last p_max iterations are plotted

x0 = 0.1; # we use the same initial value x0 for all parameters.

r = linspace(r_min, r_max, n);
pop = zeros(p_max, n);

function x =  population(r, x0, n)
# simulates n iterations of the logistic map with parameter
# r and initial value x0. The results are returned in the
# array x.
  x = zeros(n, 1);
  x(1) = x0;
  for k = 1:n-1
    x(k + 1) = r * x(k) * (1 - x(k));
  end
endfunction

for k = 1:n
  x = population(r(k), x0, t_max);
  pop(:, k) = x(t_max-p_max+1:t_max);
end

__gnuplot_set__ nokey;
__gnuplot_set__ tics out;
__gnuplot_set__ grid;
__gnuplot_set__ xlabel "{/cmmi10 r}"
__gnuplot_set__ ylabel "{/cmmi10 x}"
__gnuplot_set__ mxtics 10
__gnuplot_set__ mytics 10
__gnuplot_set__ xrange [2.8:]
__gnuplot_set__ decimalsign ","
plot(r,pop,'r.',2.995+0*(0:.1:1),0:.1:1,'b',3.545+0*(0:.1:1),0:.1:1,'b',3.45+0*(0:.1:1),0:.1:1,'b');
