This code does a simple N-body integration using arrays of doubles for storing values and a value class to make the code more readable without introducing overhead.

NBodyValClass.scala

object NBodyValClass {

	private var numBodies = 0
	private var dt = 0.0
	private var positions = Array.fill(0)(0.0)
	private var velocities = Array.fill(0)(0.0)
	private var accel = Array.fill(0)(0.0)
	private var masses = Array.fill(0)(1e-10)
	class Particle(val index: Int) extends AnyVal {
		def x = positions(index*3)
		def y = positions(index*3+1)
		def z = positions(index*3+2)
		def vx = velocities(index*3)
		def vy = velocities(index*3+1)
		def vz = velocities(index*3+2)
		def ax = accel(index*3)
		def ay = accel(index*3+1)
		def az = accel(index*3+2)
		def mass = masses(index)
		def x_=(d: Double): Unit = positions(index*3) = d
		def y_=(d: Double): Unit = positions(index*3+1) = d
		def z_=(d: Double): Unit = positions(index*3+2) = d
		def vx_=(d: Double): Unit = velocities(index*3) = d
		def vy_=(d: Double): Unit = velocities(index*3+1) = d
		def vz_=(d: Double): Unit = velocities(index*3+2) = d
		def ax_=(d: Double): Unit = accel(index*3) = d
		def ay_=(d: Double): Unit = accel(index*3+1) = d
		def az_=(d: Double): Unit = accel(index*3+2) = d
	}

	def initArrays(nb: Int, step: Double): Unit = {
		numBodies = nb
		dt = step
		positions = Array.fill(numBodies*3)(0.0)
		velocities = Array.fill(numBodies*3)(0.0)
		accel = Array.fill(numBodies*3)(0.0)
		masses = Array.fill(numBodies)(1e-10)
		masses(0) = 1.0
		for(i <-  1 until numBodies) {
			val p = new Particle(i)
			p.x = i
			p.vy = math.sqrt(1.0/i)
		}
	}

	def forSim(steps: Int): Unit = {
		require(numBodies > 0 && dt > 0.0)
		for(_ <- 1 to steps) {
			for(i <- 0 until numBodies*3) {
				accel(i) = 0.0
			}
			for {
				i <- 0 until numBodies
				j <- i+1 until numBodies
			} {
				val pi = new Particle(i)
				val pj = new Particle(j)
				val dx = pi.x-pj.x
				val dy = pi.y-pj.y
				val dz = pi.z-pj.z
				val dist = math.sqrt(dx*dx+dy*dy+dz*dz)
				val magi = pj.mass/(dist*dist*dist)
				pi.ax -= magi*dx
				pi.ay -= magi*dy
				pi.az -= magi*dz
				val magj = pi.mass/(dist*dist*dist)
				pj.ax += magj*dx
				pj.ay += magj*dy
				pj.az += magj*dz
			}
			for(i <- 0 until numBodies) {
				val p = new Particle(i)
				p.vx += p.ax*dt
				p.vy += p.ay*dt
				p.vz += p.az*dt
				p.x += p.vx*dt
				p.y += p.vy*dt
				p.z += p.vz*dt
			}
		}
	}
}