 <p>You have some acceleration function</p> <pre><code>a(p,q) where p=(x,y,z) and q=(vx,vy,vz) </code></pre> <p>Your order 1 system that can be solved via RK4 is</p> <pre><code>dotp = q dotq = a(p,q) </code></pre> <p>The stages of the RK method involve an offset of the state vector(s)</p> <pre><code>k1p = q k1q = a(p,q) p1 = p + 0.5*dt*k1p q1 = q + 0.5*dt*k1q k2p = q1 k2q = a(p1,q1) p2 = p + 0.5*dt*k2p q2 = p + 0.5*dt*k2q k3p = q2 k3q = a(p2,q2) </code></pre> <p>etc. You can either adjust the state vectors of the point <code>P</code> for each step, saving the original coordinates, or use a temporary copy of <code>P</code> to compute <code>k2, k3, k4</code>.</p>