  We present a systematic numerical study of $\theta$-dependence around
$\theta=0$ in the small-$N$ limit of $2d$ $CP^{N-1}$ models, aimed at
clarifying the possible presence of a divergent topological susceptibility in
the continuum limit. We follow a twofold strategy, based on one side on direct
simulations for $N = 2$ and $N = 3$ on lattices with correlation lengths up to
$O(10^2)$, and on the other side on the small-$N$ extrapolation of results
obtained for $N$ up to $9$. Based on that, we provide conclusive evidence for a
finite topological susceptibility at $N = 3$, with a continuum estimate $\xi^2
\chi = 0.110(5)$. On the other hand, results obtained for $N = 2$ are still
inconclusive: they are consistent with a logarithmically divergent continuum
extrapolation, but do not yet exclude a finite continuum value, $\xi^2 \chi
\sim 0.4$, with the divergence taking place for $N$ slightly below 2 in this
case. Finally, results obtained for the non-quadratic part of
$\theta$-dependence, in particular for the so-called $b_2$ coefficient, are
consistent with a $\theta$-dependence matching that of the Dilute Instanton Gas
Approximation at the point where $\xi^2 \chi$ diverges.
