BeginPackage["UVW`Lorentz`"]
(* Version 2.0, 08/15/97 *)

Lorentz::usage="
Lorentz[s,r,b] computes and plots an approximate solution of the 
Lorentz equations with parameters (s,r,b), by a Runge-Kutta method."

(* Reference : J.C. CULIOLI Introduction … Math‚matica
			    Ellipses, Paris (1991).           *)

LorentzArray::usage="
LorentzArray[matrix] plots an array of solutions of the Lorentz equations 
for the values of the parameters contained in matrix."


Begin["`Private`"]


Lorentz[s_,r_,b_,opts___]:=
      Block[{listofpoints},
	   RKuttable[f_List, y_List, y0_List, {t1_,dt_}]:=
		    Block[{k1,k2,k3,k4,yp0 = N[y0]},
			 Table[ k1 = dt N[f /. Thread[y -> yp0]];
				k2 = dt N[f /. Thread[y -> yp0 + k1/2]];
				k3 = dt N[f /. Thread[y -> yp0 + k2/2]];
				k4 = dt N[f /. Thread[y -> yp0 + k3]];
				yp0 = yp0 + (k1 + 2k2 + 2k3 + k4)/6,  
				{Round[N[t1/dt]]}]]/;
				Length[f]==Length[y]==Length[y0];

	   listofpoints= RKuttable[{- s (x-y), - x z + r x - y,x y - b z},
	      {x,y,z}, {0,1,0},{20.,0.025}];

	   Show[ Graphics3D[ {Line[listofpoints]} ], PlotRange->All, opts ]
	   ]


LorentzArray[matrix_List]:=

       Show[ GraphicsArray[
	     Table[ Lorentz[
		    matrix[[i,j,1]],matrix[[i,j,2]],matrix[[i,j,3]],
		    DisplayFunction -> Identity
			   ], 
		    {i,Length[matrix]},       
		    {j,Length[First[matrix]]}
		  ]]]



End[]

EndPackage[]