  A coupled-cluster approach for systems of $N$ bosons in external traps is
developed. In the coupled-cluster approach the exact many-body wavefunction is
obtained by applying an exponential operator $\exp{T}$ to the ground
configuration $|\phi_0>$. The natural ground configuration for bosons is, of
course, when all reside in a single orbital. Because of this simple structure
of $|\phi_0>$, the appearance of excitation operators $T=\sum_{n=1}^N T_n$ for
bosons is much simpler than for fermions. We can treat very large numbers of
bosons with coupled-cluster expansions. In a substantial part of this work, we
address the issue of size consistency for bosons and enquire whether truncated
coupled-cluster expansions are size consistent. We show that, in contrast to
the familiar situation for fermions for which coupled-cluster expansions are
size consistent, for bosons the answer to this question {\it depends} on the
choice of ground configuration. Utilizing the natural ground configuration,
working equations for the truncated coupled-cluster with $T=T_1+T_2$, i.e.,
coupled-cluster singles doubles (CCSD) are explicitly derived. Finally, an
illustrative numerical example for a condensate with up to N=10000 bosons in an
harmonic trap is provided and analyzed. The results are highly promising.
