How would you compute the $x_1,...,x_k$ that will maximise $P(x_1,...,x_k ,x_{k+1},...x_n| \mu, \Sigma)$?
if x was 2D then I think the main eigenvector of the covariance matrix at fixed $x_1$ will give you the $x_2$ that maximises $P(x_1,x_2|\mu,\Sigma)$.
However, say for n=3, you know $x_1,x_2$ and you want to know what will maximise $P(x_1,x_2,x_3 | \mu ,\Sigma)$. With $x_1,x_2$ known, you can't be on the eigenvector of the Covariance matrix anymore.
How would you do this for general n?
Since we are dealing with continuous distribution, it is meaningful to study the pdf $f(x_1,\dots,x_k|\mu,\Sigma,x_{k+1},\dots,x_n).$ Let us write the vectors $y_1=(x_1,\dots,x_k)^T$ and $y_2=(x_{k+1},\dots,x_n)^T$ For the multivariate normal,
$$\begin{pmatrix}y_1 \\y_2\end{pmatrix}\sim\mathcal{N}\left(\begin{pmatrix}\mu_1 \\\mu_2\end{pmatrix},\begin{bmatrix}\Sigma_{11} &\Sigma_{12} \\\Sigma_{21} &\Sigma_{22}\end{bmatrix}\right)$$
where $\mu_1,$ $\mu_2$ are $k\times 1$ and $(n-k)\times 1$ respectively, with $\mu=(\mu_1^T,\mu_2^T)^T.$ Similarly, $\Sigma_{11}$ is $k\times k,$ $\Sigma_{12}=\Sigma_{21}^T$ is $k\times (n-k)$ and $\Sigma_{22}$ is $(n-k)\times (n-k).$ The required conditional distribution can then be written as
$$f(y_1|\mu,\Sigma,y_2=\mathbf{a})=\mathcal{N}\left(\mu_1+\Sigma_{12}\Sigma_{22}^{-1}(\mathbf{a}-\mu_2),\Sigma_{11}-\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}\right)$$
More on this can be found here and here.
By the property of multivariate normal distribution, it is now clear that given $(x_{k+1},\dots,x_n)=\mathbf{a}^T,$ the pdf is maximized at the mean $\mu_1+\Sigma_{12}\Sigma_{22}^{-1}(\mathbf{a}-\mu_2).$