function iter_lyapunov(f,df,ri,rf,Nsteps,N,N0)
% Computation of Lyapunov exponent for a general map f = f(x,r), with given
% derivative df = df(x,r), 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: iter_lyapunov(@sinfun,@dsinfun,0.7,1.0,401,1000,200)
% (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
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 = feval(f,x,r);
        x = x1;
    end
    for i = 1:N
        x1 = feval(f,x,r);
        Lys = Lys + log(abs(feval(df,x1,r)));
        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])