function out=lyap_spec(T,trans,x0)
% Usage: lyap_spectrum(T,trans,x0)
% T is the total time and trans is the transient time
% x0 are ICs
m=length(x0);
N=T-trans; %number of iterations
% calculate orbit
% begin with initial conditions x0
x=x0;
for i=1:trans 
    [x,~]=f(x,[1.4,0.3]);
end;
% begin at this point, hopefully near attractor!
e1=zeros(m,1);
w = eye(m);
J = eye(m);
for k=1:N
    % calculate next point on trajectory and its Jacobian matrix
    [x,J]=f(x,[1.4,0.3]);
    % calculate image of unit ball under J
    % do QR decomposition to get the orthonormal vectors and
    % shrinking/stretching
    [w,R] = qr(J*w);
    % calculate stretching
    % should be e1=e1+log(norm(R(:,i)))/dt; using the R'*R matrix
    e1=e1+log(diag(R).^2)/2;
 end
% exponent is given as average
e1=e1/N; % Lyapunov exponents
out=sort(e1, 'descend');
