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

% Update v using explicit forward Euler for first iteration
n = 2;
v(:,n) = v(:,n-1) + dt*f;

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

  % Store old value of f
  f0 = f;

% Calculate rhs
  w = v(:,n-1);
  f = part_rhs(w, p);

% Update v using explicit forward Euler
  v(:,n) = v(:,n-1) + dt*(1.5*f - 0.5*f0);
  
end