PCA iteratively finds directions of greatest variance; but how to find a whole subspace with greatest variance? Principal component analysis (PCA) works like this: the first greatest variance on the first  principal component, the second greatest variance on the second principal component, and so on.
For me there is a problem with this iterative process.
What if I know I only want two principal components in order to visualize my data in 2 dimensions?
The two first PCs are not the best because the second one is the best 2nd PC but with the 1st PC they don't constitute the best couple of principal components.
Is there a way to find the "best couple" of principal components, meaning two-dimensional subspace with the greatest variance?
 A: While I completely agree with the statistical correctness of Peter Flom's answer (+1), I believe it is worth mentioning that Independent component analysis (ICA) might offer an insightful alternative. ICA provides components that are not constrained to be orthogonal with each other; this means that for some purposes ICs might be more helpful than PCs.
Looking to answer the qualitative aspect of the OP's original question: What if I know I only want two principal components in order to visualize my data in 2-dimensions ? ICs can be rather helpful as ICA tries to minimize the mutual information among the projected data.
Check for example my following shamelessly ripped-off figure (from the excellent Bayesian Reasoning and Machine Learning by David Barber (Sect 21.6) - I used this figure for a talk, I seem to have misplaced  ($\approx$ lost) my original MATLAB code, if I find it I will edit the question to append it):

Here we have a sample of two dimensional data-points (green points). The yellow lines are along the two major modes of variation in the sample defined by us during data generation. We treat them as the "true components/mode of variation".
As you see PCA's first component is indeed very close to the "true" component A. PC1 has to be the mode of maximal variation in the data. The second principal component though is constrained to be orthogonal to the first one; a condition not found in the true components. Therefore if you try to use the PCs as an alternative axis system for your data you won't get the optimal representation in terms of explanatory power (assuming that is want you want to achieve when you visualize your data; you show the data so people get the idea of what is going on). So check also ICA! It might be helpful. :)
(Once more, Peter Flom's answer is the correct answer in terms of principal components, I am huge fan of PCA don't get me wrong, just it is not always the optimal solution. As ttnphns mentioned the definition of best changes things significantly.)
(And don't be tempted to immediately equate orthogonality with statistical independence; they are related but not same; for example see the 1984 fun little paper Linearly Independent, Orthogonal, and Uncorrelated Variables by Rodger et al.)
A: The first two are the two best first two. The second one takes the first one into account. 
A: Please correct any errors you may find
Principal Components Analysis
Given data matrix $X \in \mathbb{R}^{\text{n x p}}$, where we have $n$ observations and $p$ column vectors, assume the columns of $X$ are centered with mean 0.
For any $v \in \mathbb{R}^p$, the vector $Xv \in \mathbb{R}^n$ has sample mean zero and sample variance $\frac{1}{n}(Xv)^T(Xv)$.
The first principal component direction $v_1 \in \mathbb{R}^p$ is
$v_1 = \underset{||v||_2=1}{\text{argmax}}~(Xv)^T(Xv)$
so the first principal component direction of $X$ is the unit vector $v_1 \in \mathbb{R}^p$ that maximizes the sample variance of $Xv_1 \in \mathbb{R}^p$.
The normalized first principal component score is
$\frac{(Xv_1)}{\sqrt{(Xv_1)^T(Xv_1)}}$
and the amount of variance explained by the first component is just
$\frac{\sqrt{(Xv_1)^T(Xv_1)}}{n}$
The second principal component direction is the unit vector with $v_2^Tv_1=0$, such that $Xv_2 \in \mathbb{R}^p$ has the maximum sample variance over all unit vectors orthogonal to $v_1$.
Now we can generalize further
$v_k = \underset{\underset{v^Tv_j=0, j=1,...k-1}{||v||_2=1}}{\text{argmax}}~(Xv)^T(Xv)$
for the $k^\text{th}$ principal component direction and score.
