% Maximum simulation time tmax (recommended tmax=10, tmax=20 or tmax=100)
function plotlorenz(tmax)

    % Create times with dt=0.01
    tspan=linspace(0,tmax,100*tmax);

    % Initial conditions (blue graph)
    x0=[-8 8 27];

    % Solve equation system
    [t,x]=ode45(@lorenz,tspan,x0);

    figure(1);
    clf;

    % Plot 3D trajectory
    subplot(2,2,1)
    plot3(x(:,1),x(:,2),x(:,3));

    % Plot z-value as a function of time
    subplot(2,2,2)
    plot(t,x(:,3))

    % Histogram of the distribution of the z-value
    subplot(2,2,3)
    [counts,centers]=hist(x(:,3),25);               % make histogram
    w=(centers(2)-centers(1));                      % bin width
    norm_factor = sum(counts)*w;                    % factor for normalization
    bar(centers,counts/norm_factor,1.0,'b');        % plot blue bar diagram


    % Second part: Perturb initial conditions (red graph)
    epsilon_0 = 1e-12;                              % Perturbation strength
    direction = rand(1,3)-0.5;                      % Random direction
    x02 = x0 + epsilon_0*direction/norm(direction); % Random perturbation

    % Solve equation system
    [t2,x2]=ode45(@lorenz,tspan,x02);

    % Plot 3D trajectory
    subplot(2,2,1)
    hold on;
    plot3(x2(:,1),x2(:,2),x2(:,3),'r');
    xlabel('x')
    ylabel('y')
    zlabel('z')
    title('3D Trajectory of the system in phase space')

    % Plot z-value as a function of time
    subplot(2,2,2)
    hold on;
    plot(t2,x2(:,3),'r')
    xlabel('t')
    ylabel('z(t)')
    title('z-coordinate as a function of time')
    legend('Original Solution','Perturbed Solution')

    % Histogram of the distribution of the z-value
    subplot(2,2,3)
    hold on;
    [counts,centers]=hist(x2(:,3),25);              % make histogram
    w=(centers(2)-centers(1));                      % bin width
    norm_factor = sum(counts)*w;                    % factor for normalization
    stairs(centers-0.5*w,counts/norm_factor,'r');   % plot red stair diagram
    xlabel('z')
    ylabel('p(z)')
    title('Probability distribution of the z-value')
    legend('Original Solution','Perturbed Solution')

    % Semilogarithmic plot of the error eps_z=|z_2(t)-z(t)|
    subplot(2,2,4)
    eps_z = abs((x2(:,3)-x(:,3)));                  % Absolute error in z
    lambda = 0.9;                                     % Constant lambda
    semilogy(t,eps_z,'b',t,epsilon_0*exp(lambda*t),'g')
    xlabel('t')
    ylabel('eps_z')
    title('z-Difference between the two initial conditions')
    legend('Data','f_\lambda(t)')

end

% Differential equation for the Lorenz system
function xprime = lorenz(~,x)
    sig=10;
    beta=8/3;
    rho=28;
    xprime=[-sig*x(1) + sig*x(2); rho*x(1) - x(2) - x(1)*x(3); -beta*x(3) + x(1)*x(2)];
end
