For trilinear interpolation, I think you can use Interpolations.jl:
using Interpolations
f(x,y,z) = 2x + 3y + 5z
x = y = z = -10:5.0:10
v = [f(x,y,z) for x in x, y in y, z in z]
itp = interpolate((x,y,z),v, Gridded(Linear()))
itp(1.,1.,1.) ≈ f(1.,1.,1.) # true (= 10.0)
See also interesting github discussion here.