Hessian matrix of log marginal likelihood of Gaussian Process

I'm trying to compute the exact second derivatives of log marginal likelihood of Gaussian Process for learning hyperparameters.

The log marginal likelihood and its partial derivative are given in 5 chapter of GPML(http://www.gaussianprocess.org/gpml/chapters/RW5.pdf) as

$$\log p(\mathbf{y} | \mathbf{X}, \theta) = -\dfrac{1}{2}\mathbf{y}^TK^{-1}_y \mathbf{y} -\dfrac{1}{2}\log|K_y| - \dfrac{n}{2}\log{2\pi} \tag{1}$$

$$\dfrac{\partial}{\partial \theta_j} \log p(\mathbf{y} | \mathbf{X}, \theta) = \dfrac{1}{2} \mathbf{y}^T K^{-1} \dfrac{\partial K}{\partial \theta_j} K^{-1} \mathbf{y} - \dfrac{1}{2} tr(K^{-1} \dfrac{\partial K}{\partial \theta_j}) \tag{2}$$ where

$K_y(x^p, x^q) = \sigma^2_f \exp(-\dfrac{1}{2l^2}(x_p - x_q)^T (x_p - x_q)) + \sigma^2_n\delta_{pq}$,

$\theta = (l, \sigma_f, \sigma_n)$

refferring to matrix cookbook, I'm tackling with calculating the hessian matrix of $\log p(\mathbf{y} | \mathbf{X}, \mathbf{\theta})$, but it's not easy for me to do.

Can you compute hessian matrix of $\log p(\mathbf{y} | \mathbf{X}, \mathbf{\theta})$ by hand?