
%Compute and plot the largest Lyapunov Exponent of the Lorenz System
% showplot = 1 for plot of Lyapunov exponent vs time
% showplot = 2 for additional plot of y vs x
function l = LorenzLyp(tmax,r, showplot)

% Set initial conditions and time step
tstep = 1;
x0 = [1,1,1];
v0 = [.1,.1,.1];
w0 = horzcat(x0,v0);
scalefactor = -log(norm(v0));

% Lorenz vector field parameters
sigma = 10.0;
%r = 28.0;
b = 8/3;

% Initialize data storage vectors
T=[]; W=[]; L=[];

for time = 1:tstep:tmax
    [t,w] = ode45(@LorenzOneJet,[time,time+tstep],w0,[],sigma,r,b);
    Lyp = (scalefactor + 0.5*log(sum(w(:,4:6).^2,2)))./t;
    
    % Store integration output
    T=[T; t];
    W=[W; w(:,1:3)];
    L=[L; Lyp];
    
    % Rescale deviation vector to avoid overflow and accuracy problems
    L0 = norm(w(end,4:6));
    scalefactor = scalefactor + log(L0);
    
    % Reset initial condition
    w0 = horzcat(w(end,1:3), w(end,4:6)/L0);
end  

% Create Figure and Format axes
if(showplot>=1)
    if(showplot==1)
        figure      
        % Plot Lyapunov time series
        plot(T,L);
        axis([0,tmax,min(L),1.0]);
        xlabel('t','FontSize', 14.0, 'FontName', 'Times')
        ylabel('\mu_{max}', 'Rotation',0.0, 'FontSize', 14.0, 'FontName','Times')
        %title(['Lyapunov Exponent vs Time r = ', num2str(r)])
        %grid on
    else
        figure('position',[200 200 800 300]);
        ax1=axes('position',[.07 .12 .4 .8]);
        ax2=axes('position',[.56 .12 .4 .8]);
        
        % Plot Trajectory in R3
        axes(ax1)
        plot3(W(:,1),W(:,2),W(:,3));
        xlabel('x')
        ylabel('y')
        zlabel('z')
        title('Trajectory in R^3')
        grid on
        
        % Plot Lyapunov time series
        axes(ax2);
        plot(T,L);
        axis([0,tmax,min(L),max(L)]);
        xlabel('t')
        ylabel('\mu_{max}')
        title(['Lyapunov Exponent vs Time r = ', num2str(r)])
        grid on
    end
    set(gcf,'color','white')
end

% Return last value of Lyapunov Exponent time series
l = Lyp(end);

