


%% 

% Example computation of Lyapunov exponents for the Lorenz attractor (L63 model).

%  ffun is a user supplied function that returns the differential equation
%  and the Jacobian product. With x' = f(x), and Df(x)*v the Jacobian at x, the
%  functional form is:
%     [f,Df_v] = ffun (x, v, param);
%
ffun = @fL63;  % Lorenz 63 model
x0 = [10; 10; 10];  % initial condition
param = [];   % optional constant parameters for the model (unneeded for L63)

p = 3;  % compute the p largest Lyapunov exponents

Y0 = rand(3,p);  % random initial perturbation
[Q0,R0] = mgs(Y0);

Tend = 100; %  integration time interval
tau = 0.01; %  time step
Nsteps = ceil(Tend/tau);  % number of steps

[lambda,X] = lyapunov(x0,Q0,R0,tau,Nsteps,ffun,param,Nsteps);

Lexp = cumsum(lambda,2)./repmat([1:Nsteps+1],p,1);

plot(Lexp');

Lexp(:,end)