This is essentially a "groupwise heteroskedastic" model, and it doesn't matter whether the three $Y_i$'s and the three $X_i$'s may contain different random variables (and not just different realizations of the same random variables). All that matters is that
a) The coefficient vector is assumed identical,
b) The regressor matrices have the same column dimension, and
c) The $\mathbf U_1, \mathbf U_2, \mathbf U_3$ are assumed independent.
Then it is business as usual.
To solve it, we stack our data set per equation, i.e. all observations pertaining to $Y_1$ equation first, etc, obtaining a $T\times K$ regressor matrix, $T= T_1+T_2+T_3$. We have the model
$$\mathbf y = \mathbf X\beta + \mathbf u \tag{1}$$
Mind you, we don't really care if each equation has different regressors
The model is heteroskedastic, with $T \times T$ conditional covariance matrix
$$E(\mathbf u \mathbf u' \mid \mathbf X)=\mathbf V = \text{diag}(\sigma^2_1,...,\sigma^2_2,...,\sigma^2_3,...) \tag{2}$$ with appropriate sub-dimensions. This matrix is symmetric and positive definite, so, as a general matrix-algebra result, there exists a decomposition
$$\mathbf V^{-1} = \mathbf C'\mathbf C = \text{diag}(1/\sigma^2_1,...,1/\sigma^2_2,...,1/\sigma^2_3,...) \tag{3}$$
which also gives
$$\mathbf V = \mathbf C^{-1} (\mathbf C')^{-1} \tag{4}$$
So in our case
$$\mathbf C = \text{diag}(1/\sigma_1,...,1/\sigma_2,...,1/\sigma_3,...) \tag{5}$$
where $\mathbf C$ is non-singular $T \times T$ square matrix. Consider the linear transformation of the model by premultiplication by $\mathbf C$:
$$\mathbf {\tilde y} = \mathbf C \mathbf y, \;\; \mathbf {\tilde X}= \mathbf C \mathbf X,\;\; \mathbf {\tilde u}= \mathbf C \mathbf u$$
This transformation makes the error term conditionally homoskedastic, since
$$E[\mathbf {\tilde u}\mathbf {\tilde u}'\mid \mathbf X] = \mathbf C E(\mathbf u\mathbf u'\mid \mathbf X) \mathbf C' = \mathbf C \mathbf V \mathbf C' = \mathbf C \mathbf C^{-1} (\mathbf C')^{-1} \mathbf C' = \mathbf I_T$$
and one can verify that all other properties of the linear regression model are satisfied. So running OLS on the transformed model is optimal (in the usual sense).
We will obtain
$$\hat \beta^{(tr)}_{OLS} = \left(\mathbf {\tilde X}'\mathbf {\tilde X}\right)^{-1}\mathbf {\tilde X}'\mathbf {\tilde y} = \left(\mathbf X'\mathbf C' \mathbf C\mathbf X\right)^{-1}\mathbf X' \mathbf C' \mathbf C\mathbf y = \left(\mathbf X'\mathbf V^{-1}\mathbf X\right)^{-1}\mathbf X'\mathbf V^{-1}\mathbf y $$
But this is equivalent to run GLS on the un-transformed model,
$$\beta_{GLS} = \left(\mathbf X'\mathbf V^{-1}\mathbf X\right)^{-1}\mathbf X'\mathbf V^{-1}\mathbf y $$
and its conditional variance is easy to derive as
$$\text{Var}\left[\hat \beta_{GLS}\right] = \left(\mathbf X'\mathbf V^{-1}\mathbf X\right)^{-1}$$