function [t,x,y,z] = ode_lorenz_attractor()
% Solves Lorenz attractor ODE.
% Note: do not confuse Lorenz with Lorentz (the Lorentz force, transformation, etc).
% x' = sigma*(y-x)
% y' = x*(rho - z) - y
% z' = x*y - beta*z

rho = 28; sigma = 10; beta = 8/4;

% transform system to the canonical form
%  t -> x
%  x -> y1 , sigma*(y-x)       -> f1 = sigma*(y2-y1)
%  y -> y2 , x*(rho - z) - y   -> f2 = y1*(rho - y3) -y2
%  z -> y3 , x*y - beta*z      -> f3 = y1*y2 - beta*y3

	function fvec = lorenz_attractor_forces(x,y)
		fvec(1) = sigma*(y(2)-y(1));
		fvec(2) = y(1)*(rho - y(3)) -y(2);
		fvec(3) = y(1)*y(2) - beta*y(3);

		% To make matlab happy, we need to return a column vector.
		% So, we transpose (note the dot in .')
		fvec=fvec.'; 
	end

	xspan=[0, 30];     % time interval of interest in original coordinates
	y0 = [0, 1 ,1.05]; % the initial conditions [x0,y0,z0]


	% We are using MATLAB's solver
	[x,y] = ode45(@lorenz_attractor_forces, xspan, y0); 

	% Back to the original coordinates from the ODE naming scheme
	t=x; % x is the independent variable in the ODE naming scheme
	x=y(:,1); % now x is just a coordinate
	z=y(:,3);
	y=y(:,2); % now y is just a coordinate. We need to do it the very last step.

	plot3(x,y,z);
    title('The Lorenz attractor trajectory');
	xlabel('x');
	ylabel('y');
	zlabel('z');
end
		
