Let b be your N*1 vector and f be your function values evealuated and of size N*1.
df = gradient(f,b) ; % first derivative
d2f = gradient(df,b) ; % second order gradient
The above is one method. Also you can evaluate the second derivative of f w.r.t. b and then you can substitute the values of b at the end.
Example:
syms b
f = 3*b^2+2*b+5 ;
d2f = diff(f,b,2)
