From the fact that your question stems from a finance textbook, I assume that you are dealing with Markowitz's portfolio theory. In most books, the author just focuses on the 2-asset case, but this can easily be generalized to $n$ risky assets by using matrix algebra.
Formally, let $\boldsymbol{\omega}=(\omega_1,\dots,\omega_n)^\top$ be a vector of portfolio weights and let $\boldsymbol{R}=(R_1,\dots,R_n)^\top$ be the vector of returns. I.e., $R_i$ is the return of asset $i$, which is a random variable, an $\omega_i$ is the corresponding weight, which is a scalar. The covariance matrix of $\boldsymbol{R}$ is given by:
\begin{align}
\boldsymbol{\Sigma}&=
\begin{pmatrix}
Var(R_1) & Cov(R_1,R_2) & \dots &Cov(R_1,R_n) \\
Cov(R_2,R_1) & Var(R_2) & \dots &Cov(R_2,R_n) \\
\vdots & \vdots & \ddots & \vdots \\
Cov(R_n,R_1) & Cov(R_n,R_2) & \dots & Var(R_n)
\end{pmatrix}
=\begin{pmatrix}
\sigma_1^2 & \sigma_{12} & \dots & \sigma_{1n} \\
\sigma_{21} & \sigma_2^2 & \dots &\sigma_{2n} \\
\vdots & \vdots & \ddots & \vdots \\
\sigma_{n1} & \sigma_{n2} & \dots & \sigma_n^2
\end{pmatrix}
\end{align}
Then,
$$
R_{\boldsymbol{\omega}}=\boldsymbol{\omega}^\top\boldsymbol{R}=\sum_{i=1}^n\omega_iR_i
$$
is the return of your portfolio. Basically, you want to calculate the variance of $R_{\boldsymbol{\omega}}$ and the easiest way to do this is as follows:
The expected return of the portfolio is given by:
\begin{align}
E(R_{\boldsymbol{\omega}})=E(\boldsymbol{\omega}^\top \boldsymbol{R})=\boldsymbol{\omega}^\top E(\boldsymbol{R})=\boldsymbol{\omega}^\top \boldsymbol{\mu}
\end{align}
The variance $\sigma_{\boldsymbol{\omega}}^2$ of the portfolio return is given by:
\begin{align*}
\sigma_{\boldsymbol{\omega}}^2&=E([R_{\boldsymbol{\omega}}-E(R_{\boldsymbol{\omega}})][R_{\boldsymbol{\omega}}-E(R_{\boldsymbol{\omega}})]^\top) \\
&=E([\boldsymbol{\omega}^\top\boldsymbol{R}-\boldsymbol{\omega}^\top \boldsymbol{\mu}][\boldsymbol{\omega}^\top\boldsymbol{R}-\boldsymbol{\omega}^\top \boldsymbol{\mu}]^\top) \\
&=E(\boldsymbol{\omega}^\top[\boldsymbol{R}-\boldsymbol{\mu}][\boldsymbol{\omega}^\top[\boldsymbol{R}-\boldsymbol{\mu}]]^\top) \\
&=E(\boldsymbol{\omega}^\top[\boldsymbol{R}-\boldsymbol{\mu}][\boldsymbol{R}-\boldsymbol{\mu}]^\top\boldsymbol{\omega}) \\
&=\boldsymbol{\omega}^\top E([\boldsymbol{R}-\boldsymbol{\mu}][\boldsymbol{R}-\boldsymbol{\mu}]^\top)\boldsymbol{\omega} \\
&= \boldsymbol{\omega}^\top \boldsymbol{\Sigma}\boldsymbol{\omega}\\
&=\sum_{i=1}^n\sum_{i=1}^n\omega_i\omega_jCov(R_i,R_j)
\end{align*}
Thus, to implement this in code, just specify a vector of weights, take your covariance matrix and calculate $\boldsymbol{\omega}^\top \boldsymbol{\Sigma}\boldsymbol{\omega}$.
I hope that this answers your question, if not, do not hesitate to ask further questions.