Here's a maple function for trilinear interpolation in a cube. I tried to match the notation in the wikipedia entry. I define C = (x,y,z) where x grows to the right, y grows up and z grows into the screen
tri := (x,y,z,C000,C001,C010,C011,C100,C101,C110,C111) ->
y*(z*(x*C000+(1-x)*C100)+(1-z)*(x*C010+(1-x)*C110)) +
(1-y)*(z*(x*C001+(1-x)*C101)+(1-z)*(x*C011+(1-x)*C111));