function [lambda,Xsave,Xproj] = lyapunov(x0,Q0,R0,tau,Nsteps,ffun,param,Nsave)

    x = x0;
    Q = Q0;
    R = R0;
    
    skip = floor(Nsteps/Nsave);
    
    Xsave = zeros(size(x0,1),Nsave);
    Xsave(:,1) = x0;
    
    p = size(Q,2);
    Xproj = zeros(p,Nsteps);
    
    
    lambda= zeros(size(R0,2),Nsteps+1);
    
    tol = 1e-14;

    for step=1:Nsteps,
    
        [x,Q,R] = rk4 (x,Q,R,tau,ffun,param);        
        [xdot,AQ] = ffun(x,Q,param);
        
        B = Q'*AQ;
        lambda(:,step+1) = diag(B);
     
        Xproj(:,step) = Q'*x;
        
        if mod(step,skip)==0,
            Xsave(:,step/skip+1) = x;
        end
        
%        if (norm(Q'*Q)-1)>tol,
%            disp(['reorthonogalizing in step ' num2str(step) ' err= ' num2str(norm(Q'*Q)-1) ]);
%            [Q,R] = mgs(Q*R);
        end
    end