% 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 + fa, p);
  
% Update v 
  v(:,n) = v(:,n-1) + 0.5*(fa+fb);
  
end





