If

f:= x^3 + 2yz + 5 x^2z 

is a function and if {2,0,1},{3,9,2},{8,5,6},{1,1,7} are points. I would like a Hessian of f with the four points together.

I did that with the first point like:

hessian[x_, y_, z_] = Function[{f}, D[f, {{x, y, z}}, {{x, y, z}}] ]; 
f[x_, y_, z_] :=  x^3 + 2yz + 5 x^2z;
hessian[x, y, z][f[x, y, z]]//FullSimplify
%/.{x->2,y->0,z->1} 

So, I need to apply Hessian of f at all the points at the same time. Thanks for any help