<p>I am trying to solve numerically a set of partial differential equations in three dimensions. In each of the equations the next value of the unknown in a point depends on the current value of each unknown in the closest points.</p>

<p>To write an efficient code I need to keep the points close in the three dimensions close in the (one-dimensional) memory space, so that each value is called from memory just once.</p>

<p>I was thinking of using octtrees, but I was wondering if someone knows a better method.</p>

Octtrees are the way to go.  You subdivide the array into 8 octants:
```
1 2
3 4

---

5 6
7 8
```
And then lay them out in memory in the order 1, 2, 3, 4, 5, 6, 7, 8 as above.  You repeat this recursively within each octant until you get down to some base size, probably around 128 bytes or so (this is just a guess -- make sure to profile to determine the optimal cutoff point).  This has much, much better cache coherency and locality of reference than the naive layout.