  We introduce the Nuclear Electronic All-Particle Density Matrix
Renormalization Group (NEAP-DMRG) method for solving the time-independent
Schr\"odinger equation simultaneously for electrons and other quantum species.
In contrast to already existing multicomponent approaches, in this work we
construct from the outset a multi-reference trial wave function with
stochastically optimized non-orthogonal Gaussian orbitals. By iterative
refining of the Gaussians' positions and widths, we obtain a compact
multi-reference expansion for the multicomponent wave function. We extend the
DMRG algorithm to multicomponent wave functions to take into account inter- and
intra-species correlation effects. The efficient parametrization of the total
wave function as a matrix product state allows NEAP-DMRG to accurately
approximate full configuration interaction energies of molecular systems with
more than three nuclei and twelve particles in total, which is currently a
major challenge for other multicomponent approaches. We present NEAP-DMRG
results for two few-body systems, i.e., H$_2$ and H$_3^+$, and one larger
system, namely BH$_3$
