The simplest Runge-Kutta method implementation ever (uses vector.js)

rk4.js

//Author: Anatoly Grabovsky (https://github.com/grabantot)
var euler = function(p0, v0, a_f, dt) {  //This method is here only for comparison
	var a0 = a_f(p, v)    //a=f(p0, v0)
	p0.iadd(v0.mul(dt))   //p=p0+v*dt
	v0.iadd(a0.imul(dt))  //v=p0+a*dt
}

var rk4 = function(p0, v0, a_f, dt) {
//p0 = new Vector(...), v0 = new Vector(...), a_f = function(p,v) { return new Vector(...) } }
	
	var VelAcc = function(v, a) {
		this.v = v
		this.a = a
	}
	
	var eul = function(p0, v0, a_f, dt, va) {  //Euler method for given v and a
		var p = p0.add(va.v.mul(dt))  //p=p0+v*dt
		var v = v0.add(va.a.mul(dt))  //v=p0+a*dt
		var a = a_f(p, v)   					//a=f(p, v)
		return new VelAcc(v, a)
	}

	var va0 = new VelAcc(v0, a_f(p0, v0))
	var va1 = eul(p0, v0, a_f, dt*0.5, va0)  //call Euler method 3 times
	var va2 = eul(p0, v0, a_f, dt*0.5, va1)
	var va3 = eul(p0, v0, a_f, dt*1.0, va2)
	var v = va1.v.add(va2.v).mul(2).add(va0.v).add(va3.v).div(6)  //v=(v0+2*v1+2*v2+v3)/(1+2+2+1)
	var a = va1.a.add(va2.a).mul(2).add(va0.a).add(va3.a).div(6)  //a=(a0+2*a1+2*a2+a3)/(1+2+2+1)
	p0.iadd(v.imul(dt))  //p=p0+v*dt
	v0.iadd(a.imul(dt))  //v=v0+a*dt
}