function lorenz_sim_and_traj


global SIGMA RHO BETA
SIGMA = 10.;
RHO = 28.;
BETA = 8./3.;


        xlabel('X');
        ylabel('Y');
        zlabel('Z');

        % The orbit ranges chaotically back and forth around two different points,
        % or attractors.  It is bounded, but not periodic and not convergent.
        % The numerical integration, and the display of the evolving solution,
        % are handled by the function ODE23P.

        FunFcn='lorenzeq';
        % The initial conditions below will produce good results
        %y0 = [20 5 -5];
        % Random initial conditions
        y0(1)=rand*30+5;
        y0(2)=rand*35-30;
        y0(3)=rand*40-5;
        t0=0;
        tfinal=100;
        pow = 1/3;
        tol = 0.001;

        [T,Y] = ode45(@lorenzeq,[0 tfinal],y0);
        [T2,Y2] = ode45(@lorenzeq,[0 tfinal],y0+.01*rand(3,1)');
        
        size(Y)
        size(T)
        
       figure
       set(gca,'FontSize',20)
       plot3(Y(:,1),Y(:,2),Y(:,3)) ; hold on
       plot3(Y2(:,1),Y2(:,2),Y2(:,3),'r')
       xlabel('z');ylabel('x');zlabel('y')
        %labels corresp. values x y z from class / book
       
       figure
       set(gca,'FontSize',20)
       subplot(311)
       plot(T,Y(:,1));hold on
       plot(T2,Y2(:,1),'r');hold on
        xlabel('t'); ylabel('z')
       
       subplot(312)
       plot(T,Y(:,2));hold on
              plot(T2,Y2(:,2),'r');hold on
        xlabel('t'); ylabel('x')
       
       subplot(313)
       plot(T,Y(:,3)); hold on
              plot(T2,Y2(:,3),'r');hold on
       xlabel('t'); ylabel('y')
       
       
%===============================================================================
function ydot = lorenzeq(t,y)
%LORENZEQ Equation of the Lorenz chaotic attractor.
%   ydot = lorenzeq(t,y).
%   The differential equation is written in almost linear form.

global SIGMA RHO BETA

A = [ -BETA    0     y(2)
    0  -SIGMA   SIGMA
    -y(2)   RHO    -1  ];

ydot = A*y;