You can do this using D[]. Given a function with (say) arguments x and y you can write:
D[f[x,y], {{x,y},2}]
You can use Grad[] twice, of course, but D[] will get you there directly.
Some background info here:
https://mathworld.wolfram.com/Hessian.html
On that page, they present an implementation of a Hessian operator:
HessianH[f_, x_] := D[f, {x,2}]
You'd call it like this:
HessianH[f[x,y], {x,y}]
If you don't want to have to repeat the argument list, you could do something like:
H[f_, x_]:= HessianH[Apply[f, x], x]
So now you can say:
H[f, {x,y}]
As of v12, you can access a pre-defined version of this function as follows:
ResourceFunction["HessianMatrix"][f[x,y], {x,y}]
which emits a matrix in the normal way. More info here:
https://resources.wolframcloud.com/FunctionRepository/resources/HessianMatrix