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

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

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

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