I am using smooth-particle hydrodymamics(SPH) method to model a deformable object. It is based on the procedure found in Markus Gross's point based animation book, the algorithm is attached in the image below. This algorithm uses Euler for numerical integration as seen in step 21 and 22, and I am trying to see if I can implement Runge Kutta method. However, when I look at the procedure I couldn't see how I can calculate the higher order terms. in Range Kutta y' = f(t,y) k1 = euler's method k2 = hf(t + h/2, y + k1/2) k3 = hf(t + h/2, y + k2/2) k4 = hf(t + h, y + k3) in SPH u'(t+1) = u'(t) + dt * F(t)/m I see that f(t,y) should be the acceleration, but in SPH acceleration is obtained from forces summation from neighboring elements. in this page http://gafferongames.com/game-physics/integration-basics/ it describes the f(t,y) to be mass spring approximation of acceleration a = -k * x - b*v but is this suffice to use this replacement to calculate k2, k3, and k4? what other ways can be used to approximate the acceleration?