Variance calc with weighted coefficients

I want to understand the calculation of the variance of a sum of rvs where each variable is weighted using a constant coefficient. I post the full question and answer here and I highlight the areas that I am struggling to follow.

Question.

$\displaystyle Y_1,... Y_n$ is a sample from a population with mean $\displaystyle \mu$ and variance $\displaystyle \sigma$.

The sample is not random, $\displaystyle Cov(Y_i,Y_j)={\rho}\sigma^2$ for $\displaystyle i\neq{j}$. Let $\displaystyle U=\Sigma_{i=1}^na_iY_i$.

(a) Give the condition on a constants $\displaystyle a_1,...a_n$ for U to be unbiased estimator of $\displaystyle \mu$.

(b) Under this condition, calculate MSE(U).

Answer.

(a) is not a problem, just need to calculate a if $\displaystyle E(U)=\mu$ which gives

$\displaystyle E(U)=\Sigma_{i=1}^na_iE(Y_i)=\mu\Sigma_{i=1}^na_i$ and $\displaystyle \Sigma_{i=1}^na_i=1$ is the condition.

(b) given (a) is met, then MSE(U)=Var(U):

$\displaystyle Var(U)=Var(\Sigma_{i=1}^na_iY_i)=E[\Sigma_{i=1}^na_i(Y_i-\mu)\Sigma_{j=1}^na_j(Y_j-\mu)]$ (1)

*and starting from (2) I find it tricky to follow the solution:*

$\displaystyle =\Sigma_{i=1}^n\Sigma_{j=1}^na_ia_jCov(Y_i, Y_j)$ (2)

$\displaystyle =\Sigma_{i=1}^na_i^2Var(Y_i)+\Sigma_{j=1}^n\Sigma_ {i\neq{j}}Cov(Y_i, Y_j)$ (3)

It looks to me like a sum of all entries of the variance matrix of an nx1 vector... but I cannot yet move beyond that.

Then it gets OK again, a straightforward substitution

$\displaystyle =\sigma^2\Sigma_{i=1}^na_i^2+\rho\sigma^2\Sigma_{j =1}^n\Sigma_{i\neq{j}a_ia_j}$ (4).

I'd appreciate a bit of clarification on manipulations in (2) and (3)...