% Main program: Save the program in a separate .m file and run it.
function lorenz
t=linspace(0,50,30000)'; % time variables
y0=[-1;3;4]; % Initial conditions
[t,Y] = ode45(@myLorenz,t,y0); %Invoking built-in solver 'ode45'

figure(1);
plot3(Y(:,1),Y(:,2),Y(:,3));  % Plot results
grid on;

figure(2);
plot(Y(:,1),Y(:,3));  % Plot results
grid on;

figure(3);
plot(t,Y(:,2));

function dxdt=myLorenz(t,x)
% The RHS of the Lorenz attractor
% Save this function in a separate file 'myLorenz.m'
sigma = 10;
r = 28;
b = 8/3;
dxdt=[ sigma*(x(2)-x(1));
r*x(1)-x(2)-x(1)*x(3);
x(1)*x(2)-b*x(3)];
end
end

