function lorenzRK4()  

%%%%%
% Set Lorenz
% parameters
%%%%%

a=10.0; b=8.0/3.0; c=28.0;

%%%%%%%%%%%
% Set initial condition, step length 
% and number of steps
%%%%%%%%%%%

y1(1)=0;  y2(1)=1.000; y3(1)=0.0; 

h=0.01;  

xM=40;

h2=h/2;

%%%%%%%%%
% Advance solution
%%%%%%%%%

t(1)=0;
x=t(1);

i=1;
while x<xM

i=i+1;

%%%
% Step from t(i-1) to t(i) using 
% the provisional points yt1,yt2,yt3 inbetween
%%%

% 1. Compute f(t,y)

     y1m = y1(i-1);
     y2m = y2(i-1);
     y3m = y3(i-1);
     
     
	f01 = a*(y2m-y1m);
	f02 = c*y1m - y2m - y1m*y3m;
    f03 = -b*y3m + y1m*y2m;
    
% 2. Set provisional values yt1

	yt11 = y1(i-1) + h2*f01;
	yt12 = y2(i-1) + h2*f02;
    yt13 = y3(i-1) + h2*f03;
    
% 3. Compute phase speed ft1 at provisional points

	ft11 = a*(yt12-yt11);
	ft12 = c*yt11 - yt12 - yt11*yt13;
    ft13 = -b*yt13 + yt11*yt12;
    
% 4. Set provisional values yt2

	yt21 = y1(i-1) + h2*ft11;
	yt22 = y2(i-1) + h2*ft12;
    yt23 = y3(i-1) + h2*ft13;

% 5. Compute phase speed ft2 at provisional points

	ft21 = a*(yt22-yt21);
	ft22 = c*yt21 - yt22 - yt21*yt23;
    ft23 = -b*yt23 + yt21*yt22;
    
% 6. Set provisional values yt3

	yt31 = y1(i-1) + h*ft21;
	yt32 = y2(i-1) + h*ft22;
    yt33 = y3(i-1) + h*ft23;
    
% 7. Compute phase speed ft3 at provisional points

	ft31 = a*(yt32-yt31);
	ft32 = c*yt31 - yt32 - yt31*yt33;
    ft33 = -b*yt33 + yt31*yt32;
    
% 8. Compute final phase speeds ff

	ff1 = (1/6)*( f01 + 2*ft11 + 2*ft21 + ft31 );
	ff2 = (1/6)*( f02 + 2*ft12 + 2*ft22 + ft32 );
    ff3 = (1/6)*( f03 + 2*ft13 + 2*ft23 + ft33 );
    
% 9. Set new solution values

	y1(i) = y1(i-1) + h*ff1;
	y2(i) = y2(i-1) + h*ff2;
    y3(i) = y3(i-1) + h*ff3;

% Update independent variable

	t(i) = t(i-1) + h;
        x = t(i-1);	

%x

end

%%%%%%
% Plot results
%%%%%%

% plot(t,y1,'r--');
plot(t,y1,'k');
% hold on;
% plot(t,y2,'r:');
% plot(t,y3,'b:');
xlabel('t','fontsize',20);
ylabel('x(t)','fontsize',20);
% 
% print -deps lot1.eps
% 
% figure
% 
% plot3(y1,y2,y3);
% xlabel('x'); ylabel('y'); zlabel('z')
% 
% print -deps lot2.eps
end