# Copyright 2002 Mike Perry, GPL, yada yada
#
# Got bored one night so I decided to write a Mathematica program to 
# implement a n-dimensional vector-based 4th order Runge-Kutta numerical 
# differential equation solver so's I could look at pretty spirals.
# 
# Hats off to the posterchildren, whose music stands accused of inspiring
# this monstrosity
 
NextSlope[sysv_, cpointv_] := Apply[sysv, cpointv];

VectorRungeKutta4[sysv_, initv_, steps_, stepsize_] := 
    Module[{StartPoint, Point1, Point2, Point3, Slope1, Slope2, Slope3, 
        NextPoint, t , points},
      points = {initv};
      StartPoint = initv;
      For[t = 0, t < steps, t++,
        	Slope1 = NextSlope[sysv, StartPoint];
        	Point1 = StartPoint + Slope1*stepsize/2;
        
        	Slope2 = NextSlope[sysv, Point1];
        	Point2 = StartPoint + Slope2*stepsize/2;
        
        	Slope3 = NextSlope[sysv, Point2];
        	Point3 = StartPoint + Slope3*stepsize;
        
        	Slope4 = NextSlope[sysv, Point3];
        	NextPoint = StartPoint +
            			 stepsize*(Slope1/6 + Slope2/3 + Slope3/3 + Slope4/6);
        	points = Append[points, NextPoint];
        	StartPoint = NextPoint;
        ];
      points
      ];

LinePlot3D[points_, color_, viewpoint_] := Module[{},
          
      Show[Graphics3D[{Apply[RGBColor, color], Line[points]}, {Axes -> True, 
              ViewPoint -> viewpoint, AxesLabel -> {"X", "Y", "Z"}}]];
      ];


MySigma = 10.0;
MyRho = 28.0;
MyBeta = 2.6667;
LorentzSys[x_, y_, z_] := {MySigma*(y - x), MyRho*x - y - z*x, 
      x*y - MyBeta*z};

points = VectorRungeKutta4[LorentzSys, {4, 4, 4}, 9000, .005];  
points2 = VectorRungeKutta4[LorentzSys, {4, 4, 4.01}, 9000, .005];

LinePlot3D[points, {0.5, 0, 0.5}, {1.3, -2.4, 2}];
LinePlot3D[points2, {0.5, 0, 0.5}, {1.3, -2.4, 2}];