function dy = lorenz(t,y)         % 
                                  %function to be called by 1st argument
                                  %of ode45

dy = zeros(3,1);                  %this says that dy is a vector of 3
                                  %components and initialize them to 
				  %zero

      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % The Lorenz model:                                  %
      % 				 	           % 
      %							   %
      %             dX					   %
      %            ----  = PY - PX			   %
      %             dt					   %
      %							   %
      %							   %
      %             dY  				   %
      %            ----- = -XZ + rX - Y			   %
      %             dt					   %
      %							   %
      %       						   %	
      %             dZ					   %
      %            ----- = XY - bZ  			   %
      %             dt                                     %
      %							   %
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% begin specifying set of ODEs

P = 10.0;				%Prandtl number

r = 28.0;				%ratio of Rayleigh 
					%number to 
					%critical Rayleigh
					%number

b = 8/3;				% a constant 

dy(1) = P*(y(2) - y(1));
dy(2) = -y(1)*y(3) + r*y(1) - y(2);
dy(3) = y(1)*y(2) - b*y(3);

