The $k$-th eigenfunctions are actually easy to describe: In $\mathbb{C}^{n+1}$ with unitary complex coordinates $z_0,z_1,\ldots,z_n$, write $Z = |z_0|^2+\cdots+|z_n|^2$.

Now, for a given $k\ge0$, let $H_k$ be the (real) vector space of real-valued polynomials $p(z,\bar z)$ that are homogeneous of degree $k$ in the $z$-variables and degree $k$ in the $\bar z$-variables and that are *harmonic*, i.e., they satisfy
$$
\frac{\partial^2p}{\partial z_0\partial \bar z_0}+\cdots+
\frac{\partial^2p}{\partial z_n\partial \bar z_n} = 0.
$$
One easily computes that
$$
\dim_\mathbb{R}H_k = {{n+k}\choose{k}}^2
-{{n+k-1}\choose{k-1}}^2.
$$

Then, for $p\in H_k$, the function $f_p:\mathbb{CP}^n\to\mathbb{R}$ given by
$$
f_p\bigl([z]\bigr) = \frac{p(z,\bar z)}{Z^k}
$$

is well-defined and is an eigenfunction with eigenvalue $\lambda_k$. Conversely, every eigenfunction with eigenvalue $\lambda_k$ is of this form for a unique $p\in H_k$.

**Added remark:** I was privately asked how to prove that the description I gave of the eigenvalues and eigenfunctions of $\mathbb{CP}^n$ is correct. I'm sure it's in the literature in various places, but it's easier to just give the argument, using the known description of the eigenvalues and eigenfunctions on the standard sphere. Here's the idea:

Let $S^{2n+1}\subset \mathbb{C}^{n+1}$ be the standard unit $2n{+}1$-sphere with its induced metric and let $\sigma:S^{2n+1}\to\mathbb{CP}^n$ be the quotient mapping $\sigma(z) = [z] = \mathbb{C}z$. Then $\sigma$ is a Riemannian submersion when $\mathbb{CP}^n$ is given the Fubini-Study metric (appropriately scaled; with this choice, the area of a linear $\mathbb{CP}^1\subset\mathbb{CP}^n$ is $\pi$, not $4\pi$).

If $f$ is an eigenfunction on $\mathbb{CP}^n$ with eigenvalue $\lambda_k$, then $\sigma^*(f)$ is an eigenfunction on $S^{2n+1}$ with eigenvalue $\lambda_k$, one that is, moreover, invariant under the circle action $\mathrm{e}^{i\phi}\cdot z = \mathrm{e}^{i\phi}z$ on $S^{2n+1}$.

Now, it is known that, if $\mu_d$ is the $d$-th eigenvalue of $S^{2n+1}$, then any eigenfunction $f$ with this eigenvalue is the restriction to $S^{2n+1}$ of a harmonic polynomial on $\mathbb{R}^{2n+2}$ that is homogeneous of degree $d$. We are looking for such polynomials that are also invariant under the circle action of multiplication by $\mathrm{e}^{i\phi}$. When expressed as a polynomial in the complex coordinates $z_i$ and $\bar z_i$, a polynomial that is invariant under this circle action must have the same number of $z$'s as $\bar z$'s. This can only happen if $d$ is even. Conversely, if $d=2k$, then we get the space $H_k$ as defined above. Thus, $\lambda_k = \mu_{2k}$, and the above description is justified.