The Hessian is the matrix of second partial derivatives, so yes, you should read it as $$\left (\begin {array}{cc}
\frac {\partial^2D}{\partial x^2}&\frac {\partial^2D}{\partial x\partial y}\\
\frac {\partial^2D}{\partial y\partial x}&\frac {\partial^2D}{\partial y^2}
\end {array}\right) $$How do you usually go about doing a partial differential?

Using variants of central difference methods. In the past i have approximated hessian by taking discrete derivatives twice, but that was very inaccurate. Now I am using difference of Gaussian for first derivative, but seek an efficient/best way to compute Dxx, Dxy etc given that we know first derivative was obtained via DoG, for two reasons.

A) I have a very slow computer, generating the scale space and taking DoG takes above 20 seconds for one 200x200 image (I have no gpu)
B) I need very high accuracy for this project (in terms of keypoint detection/matching), false positives and too few true positives can cause exponentially noticeable flaws.
Accuracy matters more than speed, but if a tiny increase in accuracy requires a large decrease in speed it's not worth the tradeoff.

Do it with pen and paper. You can directly populate an array with the answer to whatever precision you need.

You presumably want to convolve this with an image. And you're presumably doing that by DFTing the image, multiplying by the FT of the derivatives of the difference of Gaussian and then inverse DFTing. So you can save a lot of computation by writing down the Fourier transform of the difference of Gaussian functions and calculating that directly.

Do you know what the partial derivatives of the Gaussian function are? Can you Fourier transform them?