Notation / setting
We are considering a GP regression model:
\begin{equation}
y_i = f(x_i) + \epsilon_i
\end{equation}
where $y_i\in \mathbb{R}$,$x_i \in \mathbb{R}^d$, $f$ a Gaussian process (whose realizations are functions $f:\mathbb{R}^d\rightarrow \mathbb{R}$),
\begin{equation}
f \sim \mathrm{GP}(m(x_i), \kappa(x_i,x_j)).
\end{equation}
$n$ datapoints $(y_1,x_1), (y_2,x_2), (y_3,x_3),\ldots, (y_n,x_n)$ are given. (I use $\kappa$ to distinguish the function from the matrices $K(\cdot,\cdot)$ that contain values of $\kappa$ evaluated at certain points. The question denotes both by $K$)
How to handle $d$-dimensional inputs
The question covers computing the posterior predictive distribution for a test point (or $p$ test points) in the case $d=1$ and asks how to extend to the general $d=2,3,\ldots$.
Answer: nothing changes and the formulas from the one-dimensional case work as well in this case. Note that $m$ is then a function from $\mathbb{R}^d$ to $\mathbb{R}$ and $\kappa$ a function from $\mathbb{R}^d \times \mathbb{R}^d$ to $\mathbb{R}$.
So, for example the matrix denoted by $K(x,x)$ in the question is a $n\times n$ matrix for which $K(x,x)_{i,j} = \kappa(x_i, x_j)$ ($x_i$ and $x_j$ are $d$-dimensional but since $\kappa$ maps two $d$-dimensional vectors to a scalar, $\kappa(x_i, x_j)$ is a scalar. Similarly for $K(x,x')$ and $K(x',x')$ where $x'$ are the test points.
Thus, the dimensions of the matrices in the predictive covariance equation are $(p\times p) - (p \times n)\,(n \times n)\,(n \times p)$ independent of whether the elements of the matrices are obtained by evaluating a function $\kappa(\cdot,\cdot)$ whose arguments are $1$-dimensional of a function $\kappa(\cdot, \cdot)$ whose arguments are $d$-dimensional. In fact, the inputs could even be in some space other than $\mathbb{R}^d$ (such as if we have a categorical predictor) as long as a positive-definite covariance function can be defined.
An extra remark about the SE kernel appearing in the question
The question mentions the SE kernel
\begin{equation}
k_{f}(x_{i},x_{j}) = \sigma^{2}\exp\!\Big(-\frac{1}{2\ell^{2}}\sum_{j=1}^{q}(x_{i,j}-x_{k,j})^{2}\Big)
\end{equation}
Note that this is already a function from $\mathbb{R}^q \times \mathbb{R}^q$ to $\mathbb{R}$ (with scalar inputs there would no "$x_{i,j}$ and $x_{k,j}$" for different values of $j$. And $q$ should be $d$ if $d$ is the dimension of inputs.