Do it in pure numpy? How about copying the source of numdifftools?
What exactly is the obstacle to using numdifftools? There seem to be no licensing issues. In my experience, its a crafty piece of work; and calculating a hessian correctly, accounting for all kinds of nasty floating point issues, is no walk in the park. Even if an analytical derivative isn't too big a pain in the ass to implement, there is a good chance that what numdifftools does is more numerically stable (though in all likelihood much slower).
The only good reason for a specialized solution I can think of is speed; but be aware what you are trading it in for. If speed is your major concern though, you really cant go wrong with Theano.