function logistic_lyapunov(ri,rf,Nsteps,N,N0)
% Computation of Lyapunov exponent for the logistic map, using Nsteps
% values of r between ri and rf, starting at a fixed value of x0 (this
% could be modified) and iterating N0 steps, discarding these as
% transients, and then computing the exponent from an average over N iterates. 
% Sample usage: logistic_lyapunov(2.8,4.0,401,1000,200)
%
step = (rf-ri)/(Nsteps-1);
x = 0.3;          % Initial value
RR = []; LY = [];
r = ri-step;
for j = 1:Nsteps
    r = r + step;
    x = rand(1);
    Lys = 0;
    for i = 1:N0      % Eliminate transients - first N0 steps
        x1 = r*x*(1-x);
        x = x1;
    end
    for i = 1:N
        x1 = r*x*(1-x);
        Lys = Lys + log(abs(r*(1-2*x1)));
        x = x1;
    end
    Lyap = Lys/N;
    RR = [RR r];
    LY = [LY , Lyap];
end
figure(2)
clf;
plot([ri rf], [0 0], '-',RR',LY','b.','MarkerSize',2)
axis([ri rf -2 1])