  In this work we develop a complete variational many-body theory for a system
of $N$ trapped bosons interacting via a general two-body potential. In this
theory both the many-body basis functions {\em and} the respective expansion
coefficients are treated as variational parameters. The optimal variational
parameters are obtained {\em self-consistently} by solving a coupled system of
non-eigenvalue -- generally integro-differential -- equations to get the
one-particle functions and by diagonalizing the secular matrix problem to find
the expansion coefficients. We call this theory multi-configurational Hartree
for bosons or MCHB(M), where M specifies explicitly the number of one-particle
functions used to construct the configurations. General rules for evaluating
the matrix elements of one- and two-particle operators are derived and applied
to construct the secular Hamiltonian matrix. We discuss properties of the
derived equations. It is demonstrated that for any practical computation where
the configurational space is restricted, the description of trapped bosonic
systems strongly depends on the choice of the many-body basis set used, i.e.,
self-consistency is of great relevance. As illustrative examples we consider
bosonic systems trapped in one- and two-dimensional symmetric and asymmetric
double-well potentials. We demonstrate that self-consistency has great impact
on the predicted physical properties of the ground and excited states and show
that the lack of self-consistency may lead to physically wrong predictions. The
convergence of the general MCHB(M) scheme with a growing number M is validated
in a specific case of two bosons trapped in a symmetric double-well.
