  We investigate, by numerical simulations on a lattice, the
$\theta$-dependence of 2$d$ $CP^{N-1}$ models for a range of $N$ going from 9
to 31, combining imaginary $\theta$ and simulated tempering techniques to
improve the signal-to-noise ratio and alleviate the critical slowing down of
the topological modes. We provide continuum extrapolations for the second and
fourth order coefficients in the Taylor expansion in $\theta$ of the vacuum
energy of the theory, parameterized in terms of the topological susceptibility
$\chi$ and of the so-called $b_2$ coefficient. Those are then compared with
available analytic predictions obtained within the $1/N$ expansion, pointing
out that higher order corrections might be relevant in the explored range of
$N$, and that this fact might be related to the non-analytic behavior expected
for $N = 2$. We also consider sixth-order corrections in the $\theta$
expansion, parameterized in terms of the so-called $b_4$ coefficient: in this
case our present statistical accuracy permits to have reliable non-zero
continuum estimations only for $N \leq 11$, while for larger values we can only
set upper bounds. The sign and values obtained for $b_4$ are compared to
large-$N$ predictions, as well as to results obtained for $SU(N_c)$ Yang-Mills
theories, for which a first numerical determination is provided in this study
for the case $N_c = 2$.
