clear
load u2;

% det=16;
% b1=4;
% r=45.92;

det=10;
b1=8/3;
r=28;

V=eye(3);
S=V;b=0; 
k=20000;
for i=1:k
    J=[-det,det,0;r-u2(i,3),-1,-u2(i,1);u2(i,2),u2(i,1),-b1];
    J=eye(3)+0.002*J;
    B=J*V*S;
    [V,S,M]=svd(B);
    a_max=max(diag(S));
    S=(1/a_max)*S;
    b=b+log(a_max);
end

Lyapunov=(log(diag(S))+b)/(k*0.002)

% Lyapunov =
% 
%     0.8826
%     0.0453
%   -14.6549