function [t,x,y,z] = ode_lorenz_attractor()
% Solving Lorenz attractor ODE problem. 
% Note: do not confuse him with Lorentz (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 normal coordinates which we call xspan
	y0 = [0, 1 ,1.05] ;  % initial conditions [x0,y0,z0]


	% using matlab solver
	[x,y] = ode45(@lorenz_attractor_forces, xspan, y0); 

	% back to normal coordinates from ODE naming scheme
	t=x;        % notice that we use same name but different meaning
	x=y(:,1); 
	z=y(:,3);
	y=y(:,2); % convert ODE y matrix to normal coordinate vector

	plot3(x,y,z);
	xlabel('x');
	ylabel('y');
	zlabel('z');
end
		
