function LorenzEqs(tf)
%LorenzEqs(tf) solves the Lorenz equations to tf 
%   with x(0) = x0, y(0) = y0, z(0) = z0 given below,
%   using fourth/fifth-order Runge-Kutta or TRBDF2.
%   Here w(t) = [x(t), y(t), z(t)] & the Lorenz equations are 
%       dw/dt = F(w)
%   The Jacobian for TRBDF2 is supplied in J

tic

tspan = [0 tf];
sigma = 10; r = 28; b = 8/3;
% equilibria
eta = sqrt(b*(r-1)); % or = -sqrt(b*(r-1))
xeq = eta; 
yeq = eta; 
zeq = r-1;

w0 = [0; 1; 0]; %standard ICs
%w0 = [xeq; yeq; zeq]; % equilibrium ICs
%w0 = [0; eps; 0]; % near very unstable equilibrium at origin

% Runge-Kutta 4/5
options = odeset('RelTol',10^-3,'AbsTol',10^-6);
[t,w] = ode45(@F,tspan,w0,options);

% TRBDF2 with analytical Jacobian
%options = odeset('RelTol',10^-3,'AbsTol',10^-6,'Jacobian',@jacobian);
%[t,w] = ode23tb(@F,tspan,w0,options);

% TRBDF2 with MATLAB calculated "finite-difference" Jacobian
%options = odeset('RelTol',10^-3,'AbsTol',10^-6);
%[t,w] = ode23tb(@F,tspan,w0,options);

toc

figure
plot(t,w(:,2),'r-')
xlabel('t','FontSize',16)
ylabel('y','FontSize',16)

figure
plot(w(:,1),w(:,3),'b-',xeq,zeq,'r.',-xeq,zeq,'r.',w(1,1),w(1,3),'c.')
xlabel('x','FontSize',16)
ylabel('z','FontSize',16)
title('Lorenz attractor','FontSize',16)

figure
plot3(w(:,1),w(:,3),w(:,2),'b-',w(1,1),w(1,3),w(1,2),'c.', ...
    xeq,zeq,yeq,'r.',-xeq,zeq,-yeq,'r.')
grid on
xlabel('x','FontSize',16)
ylabel('z','FontSize',16)
zlabel('y','FontSize',16)
title('Lorenz attractor','FontSize',16)

%figure
%comet3(w(:,1),w(:,3),w(:,2))
%hold

    function wprime = F(t,w)
    wprime = [sigma*(w(2)-w(1))
         r*w(1)-w(2)-w(1)*w(3)
         w(1)*w(2)-b*w(3)];
    end

    function J = jacobian(t,w)
    J = [-sigma   sigma   0
         r-w(3)   -1      -w(1)
         w(2)     w(1)    -b   ];
    end

end