% script file Ex4.m
% DISCRETE POPULATION MODEL
% Plot the limit points of the sequence p_n
% depending on r. 
% last changes: 15.11.2001
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc;
clear all;
close all;

nr = 2000;
r  = linspace(0,4,nr);

% due to bifurcations (2-cycles) it is 
% NOT ENOUGH to plot only fixed points!

toler = 1e-8; % stopping criterion

limits   = zeros(nr,2); % due to bifurcation can have 2 limits
iter_max = 2000;

for ctr = 1:nr
  dist = 2*toler; % initialize
  p_n  = 1/3;
  p_n1 = r(ctr)*p_n*(1-p_n); % p_n+1
  iter = 0; % initialize
  while  (dist>toler)&(iter<iter_max);
    iter = iter+1;
    p_n2 = r(ctr)*p_n1*(1-p_n1); % p_n+2
  
    dist = min(abs(p_n1-p_n2),abs(p_n2-p_n));
    p_n  = p_n1; % update
    p_n1 = p_n2;
  end
  limits(ctr,1) = p_n;
  limits(ctr,2) = p_n1; % sometimes they coincide
end

figure
F = plot(r,limits(:,1),'b.',r,limits(:,2),'b.');
set(F,'markersize',0.2)    