Let's define a function for computing the Hesian of a scalar function of n variables in general. Thereafter, you can apply it to your function.
hessian[f_, vars_] := Module[{n = Length[vars]},
Table[D[D[f @@ vars, vars[[i]]], vars[[j]]], {i, 1, n}, {j, 1, n}]]
In the above function, the first input is the function whose hessian you want to compute. The second function is the list of variables of your function.
In order to test the function, we will consider a very simple function
$f(x,y) = x^2 + y^2$ . You call the function as follows
hessian[#1^2 + #2^2 &, {x, y}] // MatrixForm
or equivalently
hessian[Function[{x,y}, x^2 + y^2], {x,y}]//MatrixForm