  In this work, we show the application of the ``inverse problem'' method to
construct exact $N$ trapped soliton-like solutions of the nonlinear
Schr\"odinger or Gross-Pitaevskii equation (NLSE and GPE, respectively) in one,
two, and three spatial dimensions. This method is capable of finding the
external (confining) potentials which render specific assumed waveforms exact
solutions of the NLSE for both attractive ($g<0$) and repulsive ($g>0$)
self-interactions. For both signs of $g$, we discuss the stability with respect
to self-similar deformations and translations. For $g<0$, a critical mass
$M_c$, or equivalently the number of particles, for instabilities to arise can
often be found analytically. On the other hand, for the case with $g>0$
corresponding to repulsive self interactions which is often discussed in the
atomic physics realm of Bose-Einstein condensates (BEC), the bound solutions
are found to be always stable. For $g<0$, we also determine the critical mass
numerically by using linear stability or Bogoliubov-de Gennes analysis, and
compare these results with our analytic estimates. Various analytic forms for
the trapped $N$-soliton solutions are discussed, including sums of Gaussians or
higher-order eigenfunctions of the harmonic oscillator Hamiltonian.
