The mean-field approach to two-site Bose–Hubbard systems is well-established and leads to non-linear classical equations of motion for population imbalance and phase difference. It can, for example, be based on the representation of the solution of the time-dependent Schrödinger equation either by a single Glauber state or by a single atomic (SU(2)) coherent state [S. Wimberger et al., Phys. Rev. A 103, 023326 (2021)]. We demonstrate that quantum effects beyond the mean-field approximation are easily uncovered if, instead, a multiconfiguration ansatz with a few time-dependent SU(2) basis functions is used in the variational principle. For the case of plasma oscillations, the use of just two basis states, whose time-dependent parameters are determined variationally, already gives a good qualitative agreement of the phase space dynamics with numerically exact quantum solutions. In order to correctly account for more non-trivial effects, like macroscopic quantum self-trapping, moderately more basis states are needed. For the onset of spontaneous symmetry breaking, however, a multiplicity of 2 gives a significant improvement already. In any case, the number of variational trajectories needed for good agreement with the full quantum results is orders of magnitude smaller than that in the semi-classical case, which is based on multiple mean-field trajectories.