Gaussian Process predictive gradient with GPML

Given a Gaussian Process $$f(x) = \mathcal{GP}\big( m(x), k(x,x') \big)$$ with $$x\in\mathbb{R}^d$$, I would like to calculate $$\mathbb{E}\left[ \frac{\partial f}{\partial x_i} \right]$$; that is, the $$i^\mathrm{th}$$ coordinate of $$\mathrm{grad} f$$.

For short, I use the notation $$\partial_i \bullet$$ for the $$i^\mathrm{th}$$ gradient coordinate below. Following the GPML book by Rasmussen and Williams p.16, I think this can be done as follows: $$\mathbb{E}\left[ \partial_i f_*\ \big| X, f, X_* \right] = \partial_i \mathbb{E}\left[ f_*\ \big| X, f, X_* \right] = \partial_i K(X_*, X) K(X,X)^{-1} f$$ where $$X$$ are coordinates with known function values $$f$$, $$X_*$$ is a matrix of coordinates with unknown function values $$f_*$$, and we want the $$i^\mathrm{th}$$ coordinate of the gradient at these points $$\partial_i f_*$$.

In short, we simply need the derivative of the covariance function with respect to that coordinate. For example with the (isotropic) squared-exponential covariance: $$\partial_i k_\mathrm{SE}(x,x') = \partial_i k_\mathrm{SE}(r = \|x -x'\|_2) = -\frac{2r (x_i - x_i')}{L^2}\exp\left(-\frac{r^2}{2L^2}\right)$$

I was able to program this, but I don't want to reinvent the wheel and it seems likely that there would be a way to compute this with the GPML toolbox already. The help for covFunctions indicates that the second output dK is the "linear directional derivative function", which can apparently be called with dK(Q); but I am not sure what this means, or how to use it. Hence my questions:

Q1. Is it correct that the gradient of a GP with respect to its variables (not its parameters, as often found elsewhere in the context of GP training) is simply obtained as above, given the gradient of the covariance function?

Q2. How can one obtain these predictions using the GPML toolbox (if at all possible)?