% Initialize vector for ODE integration
v(:,1) = [u0; v0; w0; x0; y0; z0];

% Start iterative loop
for n = 2:N+1,

  % Calculate first stage 
  w0 = v(:,n-1);
  fa = dt*part_rhs(w0         , p);

  % Calculate second stage
  fb = dt*part_rhs(w0 + 0.5*fa, p);

  % Calculate third stage
  fc = dt*part_rhs(w0 + 0.5*fb, p);

  % Calculate fourth stage
  fd = dt*part_rhs(w0 +     fc, p);
  
% Update v 
  v(:,n) = v(:,n-1) + (fa + 2*fb + 2*fc + fd)/6.0;
  
end