We study a system of $A$ identical interacting bosons trapped by an external field by solving ab initio the many-body Schroedinger equation. A complete solution by using, for example, the traditional hyperspherical harmonics (HH) basis develops serious problems due to the large degeneracy of HH basis, symmetrization of the wave function, calculation of the matrix elements, etc. for large $A$. Instead of the HH basis, here we use the "potential harmonics" (PH) basis, which is a subset of HH basis. We assume that the contribution to the orbital and grand orbital [in $3(A-1)$-dimensional space of the reduced motion] quantum numbers comes only from the interacting pair. This implies inclusion of two-body correlations only and disregard of all higher-body correlations. Such an assumption is ideally suited for the Bose-Einstein condensate (BEC), which is extremely dilute. Unlike the $(3A-4)$ hyperspherical variables in HH basis, the PH basis involves only three {\it{active}} variables. It drastically reduces the number of coupled equations and calculation of the potential matrix becomes tremendously simplified, as it involves integrals over only three variables for any $A$. One can easily incorporate realistic atom-atom interactions in a straight forward manner. We study the ground and excited state properties of the condensate for both attractive and repulsive interactions for various particle number.