Edit: I have added an approach provided by @GiannisChantas. It would be great (and much appreciated) if someone could check if this approach is correct! I have also added a secondary question for which I need some help!
Let $X$, $Y$ be two multivariate random variables with mean values $\mu_X=(\mu_{X1},\ldots\mu_{Xd})^\top\in\Bbb{R}^d$, $\mu_Y=(\mu_{Y1},\ldots\mu_{Yd})^\top\in\Bbb{R}^d$ respectively. We assume that the above mean values are given as follows $$ \mu_{Xt}=\mathbf{w}_t\cdot(\mathbf{u}_x+\mathbf{e})+b_t, $$ and $$ \mu_{Yt}=\mathbf{w}_t\cdot(\mathbf{u}_y+\mathbf{e})+b_t, $$ where $\mathbf{w}_t\in\Bbb{R}^L$, $b_t\in\Bbb{R}^L$, $t=1,\ldots,d$. Moreover, $\mathbf{u}_x,\mathbf{u}_y,\mathbf{e}\in\Bbb{R}^L$. We want to evaluate the covariance matrix $$ \Sigma=\operatorname{cov}\big(X,Y\big)=E\Big[\big(X-\mu_X\big)\big(Y-\mu_Y\big)^\top\Big]. $$ If we denote the $\{\kappa,\lambda\}-th$ entry of $\Sigma$ by $\sigma_{\kappa \lambda}$, then $$ \sigma_{\kappa\lambda} = E\Big[ \big(X_\kappa-\mu_{x\kappa}\big) \big(Y_\lambda-\mu_{y\lambda}\big) \Big]. $$ After some manipulations, we have
$$ \sigma_{\kappa\lambda} = E\Big[X_\kappa Y_\lambda \Big] - E\Big[X_\kappa\mu_{y\lambda}\Big] - E\Big[Y_\lambda\mu_{x\kappa}\Big] + E\Big[\mu_{k\kappa}\mu_{y\lambda}\Big] $$
For the above individual expected values, the following hold true
$$ E\Big[X_\kappa Y_\lambda \Big] = \mu_{x\kappa}\mu_{y\lambda} + \operatorname{cov}\big(X_\kappa,Y_\lambda\big), $$
$$ E\Big[X_\kappa\mu_{y\lambda}\Big] = E\Big[X_\kappa\Big]\mu_{y\lambda} = \mu_{x\kappa}\big[\mathbf{w}_\kappa\cdot(\mathbf{u}_x+\mathbf{e})+b_\kappa\big], $$
$$ E\Big[Y_\lambda\mu_{x\kappa}\Big] = E\Big[Y_\lambda\Big]\mu_{x\kappa} = \mu_{y\lambda}\big[\mathbf{w}_\lambda\cdot(\mathbf{u}_y+\mathbf{e})+b_\lambda\big], $$
and finally,
$$ E\Big[\mu_{k\kappa}\mu_{y\lambda}\Big] = \mu_{k\kappa}\mu_{y\lambda} = \mathbf{w}_\kappa\cdot(\mathbf{u}_x+\mathbf{e}) \mathbf{w}_\lambda\cdot(\mathbf{u}_y+\mathbf{e}) + b_\lambda \mathbf{w}_\kappa\cdot(\mathbf{u}_x+\mathbf{e}) + b_\kappa \mathbf{w}_\lambda\cdot(\mathbf{u}_y+\mathbf{e}) + b_\kappa b_\lambda. $$
Consequently, the desired covariance is given as follows
$$ \sigma_{\kappa\lambda} = \mu_{x\kappa}\mu_{y\lambda} + \operatorname{cov}\big(X_\kappa,Y_\lambda\big) - \mu_{x\kappa}\big[\mathbf{w}_\lambda\cdot(\mathbf{u}_y+\mathbf{e})+b_\lambda\big] - \mu_{y\lambda}\big[\mathbf{w}_\kappa\cdot(\mathbf{u}_x+\mathbf{e})+b_\kappa\big] + \mathbf{w}_\kappa\cdot(\mathbf{u}_x+\mathbf{e}) \mathbf{w}_\lambda\cdot(\mathbf{u}_y+\mathbf{e}) + b_\lambda \mathbf{w}_\kappa\cdot(\mathbf{u}_x+\mathbf{e}) + b_\kappa \mathbf{w}_\lambda\cdot(\mathbf{u}_y+\mathbf{e}) + b_\kappa b_\lambda, $$
or
$$ \sigma_{\kappa\lambda} = \mu_{x\kappa}\mu_{y\lambda} + \operatorname{cov}\big(X_\kappa,Y_\lambda\big) + (b_\kappa-\mu_{x\kappa})\mathbf{w}_\lambda\cdot(\mathbf{u}_y+\mathbf{e}) -b_\lambda\mu_{x\kappa} + (b_\lambda-\mu_{y\lambda})\mathbf{w}_\kappa\cdot(\mathbf{u}_x+\mathbf{e}) -b_\kappa\mu_{y\lambda} + \mathbf{w}_\kappa\cdot(\mathbf{u}_x+\mathbf{e}) \mathbf{w}_\lambda\cdot(\mathbf{u}_y+\mathbf{e}) + b_\kappa b_\lambda. $$ Finally, if $X,Y$ are independent, then $\operatorname{cov}\big(X_\kappa,Y_\lambda\big)=0$, $\forall \kappa,\lambda$.
Secondary Question
Let $X$ be a random variable with mean value $\mu=(\mu_1,\ldots,\mu_t,\ldots,\mu_n)^\top$, where $\mu_t=\mathbf{w}_t\cdot(\mathbf{u}+\mathbf{e})+b_t$. What is the covariance $\sigma_{ij}=\operatorname{cov}(\mu_i,\mu_j)$? First of all, does it make sense ?
Is it $$ \sigma_{ij} = E\Big[ [\mathbf{w}_i\cdot(\mathbf{u}+\mathbf{e})+b_i] [\mathbf{w}_j\cdot(\mathbf{u}+\mathbf{e})+b_j] \Big] $$ ?
Since $\mu$ is a parameter and not a random variable, it having a covariance does not make sense. Even if you consider it as random variable, apart from the above assumptions, you have to define its probability distribution. Covariance is, like the mean value, an integral with respect to the probability distribution of the random variables. Here, you have not defined something like that.