function iter_orbit(f,ri,rf,Nsteps,N,N0)
% m-file to plot orbit diagram for general iterative map f = f(x,r), using
% Nsteps values of r between ri and rf, and plotting N iterates after
% discarding the first N0 iterates.
% Sample usage: To compute the orbit diagram for the sine function, use
% iter_orbit(@sinfun,0.7,1.0,801,400,400)
% (Warning: this could take some time; test the computer speed by using
% smaller values for Nsteps, N and N0 first...)  For long runs, the code
% can be significantly speeded up by coding the function f(x,r) into this
% m-file directly, rather than having repeated function calls using feval.
%
step = (rf-ri)/(Nsteps-1);
x = 0.3;          % Initial value - this could be modified
LT = [];  RR = [];
r = ri-step;
for j = 1:Nsteps
    r = r + step;
    x = rand(1);
    Xn = [];    
    for i = 1:N0      % Eliminate transients - first N0 steps
    x1 = feval(f,x,r);
    x = x1;
    end
    for i = 1:N
        x1 = feval(f,x,r);
        Xn = [Xn x1];
        x = x1;
    end
    RR = [RR r];
    LT = [LT , Xn'];
end
figure(2)
clf;
plot(RR',LT','b.','MarkerSize',2)
axis([ri rf 0 1])