How to calculate numerical Hessian for a multivariate Gaussian likelihood with non-linear constraints on the covariance matrix?

I'm trying to calculate standard errors for coefficient estimates in a multivariate likelihood where constraints involve covariance parameters. The Hessian output from fmincon is a scalar, and I'm having some difficulty coding it on my own.

Likelihood in multi1_llf.m, constraints in mycon.m, optimization in optim1.m