Derivative of multi-output Gaussian Process

I am working on a project where I estimate transition and measurements models for a kalman filter using Gaussian Processes. In order to linearize the models I require the Jacobian of the estimated Guassian Process.

For the single-output case this is no problem, But I am a little confused about how to do this for the multi- output case.

The posterior mean of the gaussian process would be

\begin{equation} \begin{aligned} \bar{f}_* = \mathbf{k}(\mathbf{x}_* \, \mathbf{X}) K(\mathbf{X}, \mathbf{X}) ^{-1} \mathbf{y}\\ \stackrel{\triangle}{=} \mathbf{k}(\mathbf{x}_* \, \mathbf{X}) \mathbf{\alpha} \end{aligned} \end{equation}

And the posterior Covariance

\begin{equation} \begin{aligned} C = \begin{bmatrix} k(X, X) k(X, x_*) \\ k(x_*, X) k(x_*, x_*) \end{bmatrix} \end{aligned} \end{equation}

\begin{equation} \begin{aligned} \text{covariance} = k(\mathbf{x}_*, \mathbf{x}_*) - k(\mathbf{x}_*, \mathbf{X})k(\mathbf{X}, \mathbf{X})^{-1}k(\mathbf{x}_*, \mathbf{X})^T \\ \end{aligned} \end{equation}

For the single-output case this is no problem

\begin{equation} \begin{aligned} \frac{\partial\bar{f}_*}{\partial \mathbf{x}_*} = \frac{\partial k(\mathbf{x}_* \, \mathbf{X})}{\partial \mathbf{x}_*} \mathbf{\alpha}\\ = \left[ \Lambda^{-1} \tilde{\mathbf{X}}^T_* (\mathbf{k}(\mathbf{x}_* \, \mathbf{X})^T \odot \mathbf{\alpha}) \right] \in \mathcal{R}_{D \times 1} \end{aligned} \end{equation}

However, my GP in this case maps $ f: R^2 \rightarrow R^2 $

e.g.

\begin{equation} \begin{aligned} \mathbf{x}_* = \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} \\ P\left(\begin{bmatrix} y_1\\ y_2 \end{bmatrix} \ \vert \cdots \right) = N\left( f( \begin{bmatrix} x_1\\ x_2 \end{bmatrix}) \vert \bar{f}_*, \mathbf{K} \right) \end{aligned} \end{equation}

But I am unsure how to formulate the Jacobian matrix of this function since I do not know the individual functions for each output dimension as I would in a deterministic case.

Any help would be greatly appreciated.