Lorenz system in runge kutta

Forgive me, I'm a novice when it comes to MATLAB. I'm trying to code a version of the Lorenz system of equations that produce a chaotic attractor in 3D. I'm going to have to fiddle with the parameters to produce a plot with various conditions. I've got a framework typed out but I seem to be MATLAB-illiterate. Here's my code:

    graphics_toolkit('gnuplot')
    #initial conditions:
    t(1) = 0;
    h = 0.01;
    x(1) = 2.4;
    y(1) = 1.0;
    z(1) = 0;

    a = 0.25;
    b = 4.0;
    F = 6.0;
    G = 1.0;

    #xdot = -(y^2) - (z^2) - a*x + a*F;
    #ydot = x*y - b*x*z - y + G;
    #zdot = b*x*y + x*z - z;
    U(x,y,z) = - (x^2) - (z^2) - a*x + a*F;
    V(x,y,z) = x*y - b*x*z - y + G;
    W(x,y,z) = b*x*y + x*z - z;

    for n = 1:600

      j1 = h*(U(x(n), y(n), z(n)));
      k1 = h*(V(x(n), y(n), z(n)));
      m1 = h*(W(x(n), y(n), z(n)));

      j2 = h*(U(x(n) + j1/2, y(n) + k1/2, z(n) + m1/2));
      k2 = h*(V(x(n) + j1/2, y(n) + k1/2, z(n) + m1/2));
      m2 = h*(V(x(n) + j1/2, y(n) + k1/2, z(n) + m1/2));

      j3 = h*(U(x(n) + j2/2, y(n) + k2/2, z(n) + m2/2));
      k3 = h*(V(x(n) + j2/2, y(n) + k2/2, z(n) + m2/2));
      m3 = h*(V(x(n) + j1/2, y(n) + k1/2, z(n) + m1/2));

      j4 = h*(U(x(n) + j3, y(n) + k3, z(n) + m3));
      k4 = h*(V(x(n) + j3, y(n) + k3, z(n) + m3));
      m4 = h*(W(x(n) + j3, y(n) + k3, z(n) + m3));

      t(n+1) = t(n) + h;

      x(n+1) = x(n) + (1/6) * (j1 + (2*j2) + (2*j3) + j4);

      y(n+1) = y(n) + (1/6) * (k1 + (2*k2) + (2*k3) + k4);

      z(n+1) = z(n) + (1/6) * (m1 + (2*m2) + (2*m3) + m4);

     end

    plot3(x,y,z)

It doesn't work. I'm sure I'm doing the function coding incorrectly, but I don't know any better. Help!

Also, if anyone has suggestions about the limits of the "if" statement, I'd be grateful. As I say, I'm very new to MATLAB.

Thanks in advance.

[Kermin Guo]

Check the Matlab function ODE45 with

>> type ode45.m

or

>> lookfor ode45

a simple way of testing ode set simulation without RK agorithm can be done with following code:

clear,

graphics_toolkit('gnuplot')

#initial conditions:

dt=0.01; x (1)=2.4; y(1)=1.0; z (1)=0.0;

a = 0.25; b = 4.0; F = 6.0; G = 1.0;

for i = 1 : 1000

x = x + (-y^2 - z^2 - a*x +a*F)*dt ;

y = y + (x*y -b*x*z -y +G)*dt ;

z = z + (b*x*y +x*z - z )*dt ;

xyz (i,:) = [x y z];

endfor

plot3(xyz(:,1),xyz(:,2),xyz(:,3) );